gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Insert citations into the doc (#184) - Add general citations to modules. - Add specific citations for max flow and min cost flow algorithms. - Add citations for the supported LP and MIP solvers. - Extend the main page. - Replace inproceedings entries with the journal versions. - Add a new bibtex entry about network simplex. - Remove unwanted entries.
0 6 0
default
6 files changed with 86 insertions and 93 deletions:
↑ Collapse diff ↑
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
namespace lemon {
20 20

	
21 21
/**
22 22
@defgroup datas Data Structures
23 23
This group contains the several data structures implemented in LEMON.
24 24
*/
25 25

	
26 26
/**
27 27
@defgroup graphs Graph Structures
28 28
@ingroup datas
29 29
\brief Graph structures implemented in LEMON.
30 30

	
31 31
The implementation of combinatorial algorithms heavily relies on
32 32
efficient graph implementations. LEMON offers data structures which are
33 33
planned to be easily used in an experimental phase of implementation studies,
34 34
and thereafter the program code can be made efficient by small modifications.
35 35

	
36 36
The most efficient implementation of diverse applications require the
37 37
usage of different physical graph implementations. These differences
38 38
appear in the size of graph we require to handle, memory or time usage
39 39
limitations or in the set of operations through which the graph can be
40 40
accessed.  LEMON provides several physical graph structures to meet
41 41
the diverging requirements of the possible users.  In order to save on
42 42
running time or on memory usage, some structures may fail to provide
43 43
some graph features like arc/edge or node deletion.
44 44

	
45 45
Alteration of standard containers need a very limited number of
46 46
operations, these together satisfy the everyday requirements.
47 47
In the case of graph structures, different operations are needed which do
48 48
not alter the physical graph, but gives another view. If some nodes or
49 49
arcs have to be hidden or the reverse oriented graph have to be used, then
50 50
this is the case. It also may happen that in a flow implementation
51 51
the residual graph can be accessed by another algorithm, or a node-set
52 52
is to be shrunk for another algorithm.
53 53
LEMON also provides a variety of graphs for these requirements called
54 54
\ref graph_adaptors "graph adaptors". Adaptors cannot be used alone but only
55 55
in conjunction with other graph representations.
56 56

	
57 57
You are free to use the graph structure that fit your requirements
58 58
the best, most graph algorithms and auxiliary data structures can be used
59 59
with any graph structure.
60 60

	
61 61
<b>See also:</b> \ref graph_concepts "Graph Structure Concepts".
62 62
*/
63 63

	
64 64
/**
65 65
@defgroup graph_adaptors Adaptor Classes for Graphs
66 66
@ingroup graphs
67 67
\brief Adaptor classes for digraphs and graphs
68 68

	
69 69
This group contains several useful adaptor classes for digraphs and graphs.
70 70

	
71 71
The main parts of LEMON are the different graph structures, generic
72 72
graph algorithms, graph concepts, which couple them, and graph
73 73
adaptors. While the previous notions are more or less clear, the
74 74
latter one needs further explanation. Graph adaptors are graph classes
75 75
which serve for considering graph structures in different ways.
76 76

	
77 77
A short example makes this much clearer.  Suppose that we have an
78 78
instance \c g of a directed graph type, say ListDigraph and an algorithm
79 79
\code
80 80
template <typename Digraph>
81 81
int algorithm(const Digraph&);
82 82
\endcode
83 83
is needed to run on the reverse oriented graph.  It may be expensive
84 84
(in time or in memory usage) to copy \c g with the reversed
85 85
arcs.  In this case, an adaptor class is used, which (according
86 86
to LEMON \ref concepts::Digraph "digraph concepts") works as a digraph.
87 87
The adaptor uses the original digraph structure and digraph operations when
88 88
methods of the reversed oriented graph are called.  This means that the adaptor
89 89
have minor memory usage, and do not perform sophisticated algorithmic
90 90
actions.  The purpose of it is to give a tool for the cases when a
91 91
graph have to be used in a specific alteration.  If this alteration is
92 92
obtained by a usual construction like filtering the node or the arc set or
93 93
considering a new orientation, then an adaptor is worthwhile to use.
94 94
To come back to the reverse oriented graph, in this situation
95 95
\code
96 96
template<typename Digraph> class ReverseDigraph;
97 97
\endcode
98 98
template class can be used. The code looks as follows
99 99
\code
100 100
ListDigraph g;
101 101
ReverseDigraph<ListDigraph> rg(g);
102 102
int result = algorithm(rg);
103 103
\endcode
104 104
During running the algorithm, the original digraph \c g is untouched.
105 105
This techniques give rise to an elegant code, and based on stable
106 106
graph adaptors, complex algorithms can be implemented easily.
107 107

	
108 108
In flow, circulation and matching problems, the residual
109 109
graph is of particular importance. Combining an adaptor implementing
110 110
this with shortest path algorithms or minimum mean cycle algorithms,
111 111
a range of weighted and cardinality optimization algorithms can be
112 112
obtained. For other examples, the interested user is referred to the
113 113
detailed documentation of particular adaptors.
114 114

	
115 115
The behavior of graph adaptors can be very different. Some of them keep
116 116
capabilities of the original graph while in other cases this would be
117 117
meaningless. This means that the concepts that they meet depend
118 118
on the graph adaptor, and the wrapped graph.
119 119
For example, if an arc of a reversed digraph is deleted, this is carried
120 120
out by deleting the corresponding arc of the original digraph, thus the
121 121
adaptor modifies the original digraph.
122 122
However in case of a residual digraph, this operation has no sense.
123 123

	
124 124
Let us stand one more example here to simplify your work.
125 125
ReverseDigraph has constructor
126 126
\code
127 127
ReverseDigraph(Digraph& digraph);
128 128
\endcode
129 129
This means that in a situation, when a <tt>const %ListDigraph&</tt>
130 130
reference to a graph is given, then it have to be instantiated with
131 131
<tt>Digraph=const %ListDigraph</tt>.
132 132
\code
133 133
int algorithm1(const ListDigraph& g) {
134 134
  ReverseDigraph<const ListDigraph> rg(g);
135 135
  return algorithm2(rg);
136 136
}
137 137
\endcode
138 138
*/
139 139

	
140 140
/**
141 141
@defgroup maps Maps
142 142
@ingroup datas
143 143
\brief Map structures implemented in LEMON.
144 144

	
145 145
This group contains the map structures implemented in LEMON.
146 146

	
147 147
LEMON provides several special purpose maps and map adaptors that e.g. combine
148 148
new maps from existing ones.
149 149

	
150 150
<b>See also:</b> \ref map_concepts "Map Concepts".
151 151
*/
152 152

	
153 153
/**
154 154
@defgroup graph_maps Graph Maps
155 155
@ingroup maps
156 156
\brief Special graph-related maps.
157 157

	
158 158
This group contains maps that are specifically designed to assign
159 159
values to the nodes and arcs/edges of graphs.
160 160

	
161 161
If you are looking for the standard graph maps (\c NodeMap, \c ArcMap,
162 162
\c EdgeMap), see the \ref graph_concepts "Graph Structure Concepts".
163 163
*/
164 164

	
165 165
/**
166 166
\defgroup map_adaptors Map Adaptors
167 167
\ingroup maps
168 168
\brief Tools to create new maps from existing ones
169 169

	
170 170
This group contains map adaptors that are used to create "implicit"
171 171
maps from other maps.
172 172

	
173 173
Most of them are \ref concepts::ReadMap "read-only maps".
174 174
They can make arithmetic and logical operations between one or two maps
175 175
(negation, shifting, addition, multiplication, logical 'and', 'or',
176 176
'not' etc.) or e.g. convert a map to another one of different Value type.
177 177

	
178 178
The typical usage of this classes is passing implicit maps to
179 179
algorithms.  If a function type algorithm is called then the function
180 180
type map adaptors can be used comfortable. For example let's see the
181 181
usage of map adaptors with the \c graphToEps() function.
182 182
\code
183 183
  Color nodeColor(int deg) {
184 184
    if (deg >= 2) {
185 185
      return Color(0.5, 0.0, 0.5);
186 186
    } else if (deg == 1) {
187 187
      return Color(1.0, 0.5, 1.0);
188 188
    } else {
189 189
      return Color(0.0, 0.0, 0.0);
190 190
    }
191 191
  }
192 192

	
193 193
  Digraph::NodeMap<int> degree_map(graph);
194 194

	
195 195
  graphToEps(graph, "graph.eps")
196 196
    .coords(coords).scaleToA4().undirected()
197 197
    .nodeColors(composeMap(functorToMap(nodeColor), degree_map))
198 198
    .run();
199 199
\endcode
200 200
The \c functorToMap() function makes an \c int to \c Color map from the
201 201
\c nodeColor() function. The \c composeMap() compose the \c degree_map
202 202
and the previously created map. The composed map is a proper function to
203 203
get the color of each node.
204 204

	
205 205
The usage with class type algorithms is little bit harder. In this
206 206
case the function type map adaptors can not be used, because the
207 207
function map adaptors give back temporary objects.
208 208
\code
209 209
  Digraph graph;
210 210

	
211 211
  typedef Digraph::ArcMap<double> DoubleArcMap;
212 212
  DoubleArcMap length(graph);
213 213
  DoubleArcMap speed(graph);
214 214

	
215 215
  typedef DivMap<DoubleArcMap, DoubleArcMap> TimeMap;
216 216
  TimeMap time(length, speed);
217 217

	
218 218
  Dijkstra<Digraph, TimeMap> dijkstra(graph, time);
219 219
  dijkstra.run(source, target);
220 220
\endcode
221 221
We have a length map and a maximum speed map on the arcs of a digraph.
222 222
The minimum time to pass the arc can be calculated as the division of
223 223
the two maps which can be done implicitly with the \c DivMap template
224 224
class. We use the implicit minimum time map as the length map of the
225 225
\c Dijkstra algorithm.
226 226
*/
227 227

	
228 228
/**
229 229
@defgroup paths Path Structures
230 230
@ingroup datas
231 231
\brief %Path structures implemented in LEMON.
232 232

	
233 233
This group contains the path structures implemented in LEMON.
234 234

	
235 235
LEMON provides flexible data structures to work with paths.
236 236
All of them have similar interfaces and they can be copied easily with
237 237
assignment operators and copy constructors. This makes it easy and
238 238
efficient to have e.g. the Dijkstra algorithm to store its result in
239 239
any kind of path structure.
240 240

	
241 241
\sa \ref concepts::Path "Path concept"
242 242
*/
243 243

	
244 244
/**
245 245
@defgroup heaps Heap Structures
246 246
@ingroup datas
247 247
\brief %Heap structures implemented in LEMON.
248 248

	
249 249
This group contains the heap structures implemented in LEMON.
250 250

	
251 251
LEMON provides several heap classes. They are efficient implementations
252 252
of the abstract data type \e priority \e queue. They store items with
253 253
specified values called \e priorities in such a way that finding and
254 254
removing the item with minimum priority are efficient.
255 255
The basic operations are adding and erasing items, changing the priority
256 256
of an item, etc.
257 257

	
258 258
Heaps are crucial in several algorithms, such as Dijkstra and Prim.
259 259
The heap implementations have the same interface, thus any of them can be
260 260
used easily in such algorithms.
261 261

	
262 262
\sa \ref concepts::Heap "Heap concept"
263 263
*/
264 264

	
265 265
/**
266 266
@defgroup matrices Matrices
267 267
@ingroup datas
268 268
\brief Two dimensional data storages implemented in LEMON.
269 269

	
270 270
This group contains two dimensional data storages implemented in LEMON.
271 271
*/
272 272

	
273 273
/**
274 274
@defgroup auxdat Auxiliary Data Structures
275 275
@ingroup datas
276 276
\brief Auxiliary data structures implemented in LEMON.
277 277

	
278 278
This group contains some data structures implemented in LEMON in
279 279
order to make it easier to implement combinatorial algorithms.
280 280
*/
281 281

	
282 282
/**
283 283
@defgroup geomdat Geometric Data Structures
284 284
@ingroup auxdat
285 285
\brief Geometric data structures implemented in LEMON.
286 286

	
287 287
This group contains geometric data structures implemented in LEMON.
288 288

	
289 289
 - \ref lemon::dim2::Point "dim2::Point" implements a two dimensional
290 290
   vector with the usual operations.
291 291
 - \ref lemon::dim2::Box "dim2::Box" can be used to determine the
292 292
   rectangular bounding box of a set of \ref lemon::dim2::Point
293 293
   "dim2::Point"'s.
294 294
*/
295 295

	
296 296
/**
297 297
@defgroup matrices Matrices
298 298
@ingroup auxdat
299 299
\brief Two dimensional data storages implemented in LEMON.
300 300

	
301 301
This group contains two dimensional data storages implemented in LEMON.
302 302
*/
303 303

	
304 304
/**
305 305
@defgroup algs Algorithms
306 306
\brief This group contains the several algorithms
307 307
implemented in LEMON.
308 308

	
309 309
This group contains the several algorithms
310 310
implemented in LEMON.
311 311
*/
312 312

	
313 313
/**
314 314
@defgroup search Graph Search
315 315
@ingroup algs
316 316
\brief Common graph search algorithms.
317 317

	
318 318
This group contains the common graph search algorithms, namely
319
\e breadth-first \e search (BFS) and \e depth-first \e search (DFS).
319
\e breadth-first \e search (BFS) and \e depth-first \e search (DFS)
320
\ref clrs01algorithms.
320 321
*/
321 322

	
322 323
/**
323 324
@defgroup shortest_path Shortest Path Algorithms
324 325
@ingroup algs
325 326
\brief Algorithms for finding shortest paths.
326 327

	
327
This group contains the algorithms for finding shortest paths in digraphs.
328
This group contains the algorithms for finding shortest paths in digraphs
329
\ref clrs01algorithms.
328 330

	
329 331
 - \ref Dijkstra algorithm for finding shortest paths from a source node
330 332
   when all arc lengths are non-negative.
331 333
 - \ref BellmanFord "Bellman-Ford" algorithm for finding shortest paths
332 334
   from a source node when arc lenghts can be either positive or negative,
333 335
   but the digraph should not contain directed cycles with negative total
334 336
   length.
335 337
 - \ref FloydWarshall "Floyd-Warshall" and \ref Johnson "Johnson" algorithms
336 338
   for solving the \e all-pairs \e shortest \e paths \e problem when arc
337 339
   lenghts can be either positive or negative, but the digraph should
338 340
   not contain directed cycles with negative total length.
339 341
 - \ref Suurballe A successive shortest path algorithm for finding
340 342
   arc-disjoint paths between two nodes having minimum total length.
341 343
*/
342 344

	
343 345
/**
344 346
@defgroup spantree Minimum Spanning Tree Algorithms
345 347
@ingroup algs
346 348
\brief Algorithms for finding minimum cost spanning trees and arborescences.
347 349

	
348 350
This group contains the algorithms for finding minimum cost spanning
349
trees and arborescences.
351
trees and arborescences \ref clrs01algorithms.
350 352
*/
351 353

	
352 354
/**
353 355
@defgroup max_flow Maximum Flow Algorithms
354 356
@ingroup algs
355 357
\brief Algorithms for finding maximum flows.
356 358

	
357 359
This group contains the algorithms for finding maximum flows and
358
feasible circulations.
360
feasible circulations \ref clrs01algorithms, \ref amo93networkflows.
359 361

	
360 362
The \e maximum \e flow \e problem is to find a flow of maximum value between
361 363
a single source and a single target. Formally, there is a \f$G=(V,A)\f$
362 364
digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and
363 365
\f$s, t \in V\f$ source and target nodes.
364 366
A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the
365 367
following optimization problem.
366 368

	
367 369
\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f]
368 370
\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu)
369 371
    \quad \forall u\in V\setminus\{s,t\} \f]
370 372
\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f]
371 373

	
372 374
LEMON contains several algorithms for solving maximum flow problems:
373
- \ref EdmondsKarp Edmonds-Karp algorithm.
374
- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm.
375
- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees.
376
- \ref GoldbergTarjan Preflow push-relabel algorithm with dynamic trees.
375
- \ref EdmondsKarp Edmonds-Karp algorithm
376
  \ref edmondskarp72theoretical.
377
- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm
378
  \ref goldberg88newapproach.
379
- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees
380
  \ref dinic70algorithm, \ref sleator83dynamic.
381
- \ref GoldbergTarjan !Preflow push-relabel algorithm with dynamic trees
382
  \ref goldberg88newapproach, \ref sleator83dynamic.
377 383

	
378
In most cases the \ref Preflow "Preflow" algorithm provides the
384
In most cases the \ref Preflow algorithm provides the
379 385
fastest method for computing a maximum flow. All implementations
380 386
also provide functions to query the minimum cut, which is the dual
381 387
problem of maximum flow.
382 388

	
383 389
\ref Circulation is a preflow push-relabel algorithm implemented directly 
384 390
for finding feasible circulations, which is a somewhat different problem,
385 391
but it is strongly related to maximum flow.
386 392
For more information, see \ref Circulation.
387 393
*/
388 394

	
389 395
/**
390 396
@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms
391 397
@ingroup algs
392 398

	
393 399
\brief Algorithms for finding minimum cost flows and circulations.
394 400

	
395 401
This group contains the algorithms for finding minimum cost flows and
396
circulations. For more information about this problem and its dual
397
solution see \ref min_cost_flow "Minimum Cost Flow Problem".
402
circulations \ref amo93networkflows. For more information about this
403
problem and its dual solution, see \ref min_cost_flow
404
"Minimum Cost Flow Problem".
398 405

	
399 406
LEMON contains several algorithms for this problem.
400 407
 - \ref NetworkSimplex Primal Network Simplex algorithm with various
401
   pivot strategies.
408
   pivot strategies \ref dantzig63linearprog, \ref kellyoneill91netsimplex.
402 409
 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
403
   cost scaling.
410
   cost scaling \ref goldberg90approximation, \ref goldberg97efficient,
411
   \ref bunnagel98efficient.
404 412
 - \ref CapacityScaling Successive Shortest %Path algorithm with optional
405
   capacity scaling.
406
 - \ref CancelAndTighten The Cancel and Tighten algorithm.
407
 - \ref CycleCanceling Cycle-Canceling algorithms.
413
   capacity scaling \ref edmondskarp72theoretical.
414
 - \ref CancelAndTighten The Cancel and Tighten algorithm
415
   \ref goldberg89cyclecanceling.
416
 - \ref CycleCanceling Cycle-Canceling algorithms
417
   \ref klein67primal, \ref goldberg89cyclecanceling.
408 418

	
409 419
In general NetworkSimplex is the most efficient implementation,
410 420
but in special cases other algorithms could be faster.
411 421
For example, if the total supply and/or capacities are rather small,
412 422
CapacityScaling is usually the fastest algorithm (without effective scaling).
413 423
*/
414 424

	
415 425
/**
416 426
@defgroup min_cut Minimum Cut Algorithms
417 427
@ingroup algs
418 428

	
419 429
\brief Algorithms for finding minimum cut in graphs.
420 430

	
421 431
This group contains the algorithms for finding minimum cut in graphs.
422 432

	
423 433
The \e minimum \e cut \e problem is to find a non-empty and non-complete
424 434
\f$X\f$ subset of the nodes with minimum overall capacity on
425 435
outgoing arcs. Formally, there is a \f$G=(V,A)\f$ digraph, a
426 436
\f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function. The minimum
427 437
cut is the \f$X\f$ solution of the next optimization problem:
428 438

	
429 439
\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}}
430 440
    \sum_{uv\in A: u\in X, v\not\in X}cap(uv) \f]
431 441

	
432 442
LEMON contains several algorithms related to minimum cut problems:
433 443

	
434 444
- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut
435 445
  in directed graphs.
436 446
- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for
437 447
  calculating minimum cut in undirected graphs.
438 448
- \ref GomoryHu "Gomory-Hu tree computation" for calculating
439 449
  all-pairs minimum cut in undirected graphs.
440 450

	
441 451
If you want to find minimum cut just between two distinict nodes,
442 452
see the \ref max_flow "maximum flow problem".
443 453
*/
444 454

	
445 455
/**
446 456
@defgroup matching Matching Algorithms
447 457
@ingroup algs
448 458
\brief Algorithms for finding matchings in graphs and bipartite graphs.
449 459

	
450 460
This group contains the algorithms for calculating
451 461
matchings in graphs and bipartite graphs. The general matching problem is
452 462
finding a subset of the edges for which each node has at most one incident
453 463
edge.
454 464

	
455 465
There are several different algorithms for calculate matchings in
456 466
graphs.  The matching problems in bipartite graphs are generally
457 467
easier than in general graphs. The goal of the matching optimization
458 468
can be finding maximum cardinality, maximum weight or minimum cost
459 469
matching. The search can be constrained to find perfect or
460 470
maximum cardinality matching.
461 471

	
462 472
The matching algorithms implemented in LEMON:
463 473
- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm
464 474
  for calculating maximum cardinality matching in bipartite graphs.
465 475
- \ref PrBipartiteMatching Push-relabel algorithm
466 476
  for calculating maximum cardinality matching in bipartite graphs.
467 477
- \ref MaxWeightedBipartiteMatching
468 478
  Successive shortest path algorithm for calculating maximum weighted
469 479
  matching and maximum weighted bipartite matching in bipartite graphs.
470 480
- \ref MinCostMaxBipartiteMatching
471 481
  Successive shortest path algorithm for calculating minimum cost maximum
472 482
  matching in bipartite graphs.
473 483
- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating
474 484
  maximum cardinality matching in general graphs.
475 485
- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating
476 486
  maximum weighted matching in general graphs.
477 487
- \ref MaxWeightedPerfectMatching
478 488
  Edmond's blossom shrinking algorithm for calculating maximum weighted
479 489
  perfect matching in general graphs.
480 490

	
481 491
\image html bipartite_matching.png
482 492
\image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth
483 493
*/
484 494

	
485 495
/**
486 496
@defgroup graph_properties Connectivity and Other Graph Properties
487 497
@ingroup algs
488 498
\brief Algorithms for discovering the graph properties
489 499

	
490 500
This group contains the algorithms for discovering the graph properties
491 501
like connectivity, bipartiteness, euler property, simplicity etc.
492 502

	
493 503
\image html connected_components.png
494 504
\image latex connected_components.eps "Connected components" width=\textwidth
495 505
*/
496 506

	
497 507
/**
498 508
@defgroup planar Planarity Embedding and Drawing
499 509
@ingroup algs
500 510
\brief Algorithms for planarity checking, embedding and drawing
501 511

	
502 512
This group contains the algorithms for planarity checking,
503 513
embedding and drawing.
504 514

	
505 515
\image html planar.png
506 516
\image latex planar.eps "Plane graph" width=\textwidth
507 517
*/
508 518

	
509 519
/**
510 520
@defgroup approx Approximation Algorithms
511 521
@ingroup algs
512 522
\brief Approximation algorithms.
513 523

	
514 524
This group contains the approximation and heuristic algorithms
515 525
implemented in LEMON.
516 526
*/
517 527

	
518 528
/**
519 529
@defgroup auxalg Auxiliary Algorithms
520 530
@ingroup algs
521 531
\brief Auxiliary algorithms implemented in LEMON.
522 532

	
523 533
This group contains some algorithms implemented in LEMON
524 534
in order to make it easier to implement complex algorithms.
525 535
*/
526 536

	
527 537
/**
528 538
@defgroup gen_opt_group General Optimization Tools
529 539
\brief This group contains some general optimization frameworks
530 540
implemented in LEMON.
531 541

	
532 542
This group contains some general optimization frameworks
533 543
implemented in LEMON.
534 544
*/
535 545

	
536 546
/**
537
@defgroup lp_group Lp and Mip Solvers
547
@defgroup lp_group LP and MIP Solvers
538 548
@ingroup gen_opt_group
539
\brief Lp and Mip solver interfaces for LEMON.
549
\brief LP and MIP solver interfaces for LEMON.
540 550

	
541
This group contains Lp and Mip solver interfaces for LEMON. The
542
various LP solvers could be used in the same manner with this
543
interface.
551
This group contains LP and MIP solver interfaces for LEMON.
552
Various LP solvers could be used in the same manner with this
553
high-level interface.
554

	
555
The currently supported solvers are \ref glpk, \ref clp, \ref cbc,
556
\ref cplex, \ref soplex.
544 557
*/
545 558

	
546 559
/**
547 560
@defgroup lp_utils Tools for Lp and Mip Solvers
548 561
@ingroup lp_group
549 562
\brief Helper tools to the Lp and Mip solvers.
550 563

	
551 564
This group adds some helper tools to general optimization framework
552 565
implemented in LEMON.
553 566
*/
554 567

	
555 568
/**
556 569
@defgroup metah Metaheuristics
557 570
@ingroup gen_opt_group
558 571
\brief Metaheuristics for LEMON library.
559 572

	
560 573
This group contains some metaheuristic optimization tools.
561 574
*/
562 575

	
563 576
/**
564 577
@defgroup utils Tools and Utilities
565 578
\brief Tools and utilities for programming in LEMON
566 579

	
567 580
Tools and utilities for programming in LEMON.
568 581
*/
569 582

	
570 583
/**
571 584
@defgroup gutils Basic Graph Utilities
572 585
@ingroup utils
573 586
\brief Simple basic graph utilities.
574 587

	
575 588
This group contains some simple basic graph utilities.
576 589
*/
577 590

	
578 591
/**
579 592
@defgroup misc Miscellaneous Tools
580 593
@ingroup utils
581 594
\brief Tools for development, debugging and testing.
582 595

	
583 596
This group contains several useful tools for development,
584 597
debugging and testing.
585 598
*/
586 599

	
587 600
/**
588 601
@defgroup timecount Time Measuring and Counting
589 602
@ingroup misc
590 603
\brief Simple tools for measuring the performance of algorithms.
591 604

	
592 605
This group contains simple tools for measuring the performance
593 606
of algorithms.
594 607
*/
595 608

	
596 609
/**
597 610
@defgroup exceptions Exceptions
598 611
@ingroup utils
599 612
\brief Exceptions defined in LEMON.
600 613

	
601 614
This group contains the exceptions defined in LEMON.
602 615
*/
603 616

	
604 617
/**
605 618
@defgroup io_group Input-Output
606 619
\brief Graph Input-Output methods
607 620

	
608 621
This group contains the tools for importing and exporting graphs
609 622
and graph related data. Now it supports the \ref lgf-format
610 623
"LEMON Graph Format", the \c DIMACS format and the encapsulated
611 624
postscript (EPS) format.
612 625
*/
613 626

	
614 627
/**
615 628
@defgroup lemon_io LEMON Graph Format
616 629
@ingroup io_group
617 630
\brief Reading and writing LEMON Graph Format.
618 631

	
619 632
This group contains methods for reading and writing
620 633
\ref lgf-format "LEMON Graph Format".
621 634
*/
622 635

	
623 636
/**
624 637
@defgroup eps_io Postscript Exporting
625 638
@ingroup io_group
626 639
\brief General \c EPS drawer and graph exporter
627 640

	
628 641
This group contains general \c EPS drawing methods and special
629 642
graph exporting tools.
630 643
*/
631 644

	
632 645
/**
633 646
@defgroup dimacs_group DIMACS Format
634 647
@ingroup io_group
635 648
\brief Read and write files in DIMACS format
636 649

	
637 650
Tools to read a digraph from or write it to a file in DIMACS format data.
638 651
*/
639 652

	
640 653
/**
641 654
@defgroup nauty_group NAUTY Format
642 655
@ingroup io_group
643 656
\brief Read \e Nauty format
644 657

	
645 658
Tool to read graphs from \e Nauty format data.
646 659
*/
647 660

	
648 661
/**
649 662
@defgroup concept Concepts
650 663
\brief Skeleton classes and concept checking classes
651 664

	
652 665
This group contains the data/algorithm skeletons and concept checking
653 666
classes implemented in LEMON.
654 667

	
655 668
The purpose of the classes in this group is fourfold.
656 669

	
657 670
- These classes contain the documentations of the %concepts. In order
658 671
  to avoid document multiplications, an implementation of a concept
659 672
  simply refers to the corresponding concept class.
660 673

	
661 674
- These classes declare every functions, <tt>typedef</tt>s etc. an
662 675
  implementation of the %concepts should provide, however completely
663 676
  without implementations and real data structures behind the
664 677
  interface. On the other hand they should provide nothing else. All
665 678
  the algorithms working on a data structure meeting a certain concept
666 679
  should compile with these classes. (Though it will not run properly,
667 680
  of course.) In this way it is easily to check if an algorithm
668 681
  doesn't use any extra feature of a certain implementation.
669 682

	
670 683
- The concept descriptor classes also provide a <em>checker class</em>
671 684
  that makes it possible to check whether a certain implementation of a
672 685
  concept indeed provides all the required features.
673 686

	
674 687
- Finally, They can serve as a skeleton of a new implementation of a concept.
675 688
*/
676 689

	
677 690
/**
678 691
@defgroup graph_concepts Graph Structure Concepts
679 692
@ingroup concept
680 693
\brief Skeleton and concept checking classes for graph structures
681 694

	
682 695
This group contains the skeletons and concept checking classes of
683 696
graph structures.
684 697
*/
685 698

	
686 699
/**
687 700
@defgroup map_concepts Map Concepts
688 701
@ingroup concept
689 702
\brief Skeleton and concept checking classes for maps
690 703

	
691 704
This group contains the skeletons and concept checking classes of maps.
692 705
*/
693 706

	
694 707
/**
695 708
@defgroup tools Standalone Utility Applications
696 709

	
697 710
Some utility applications are listed here.
698 711

	
699 712
The standard compilation procedure (<tt>./configure;make</tt>) will compile
700 713
them, as well.
701 714
*/
702 715

	
703 716
/**
704 717
\anchor demoprograms
705 718

	
706 719
@defgroup demos Demo Programs
707 720

	
708 721
Some demo programs are listed here. Their full source codes can be found in
709 722
the \c demo subdirectory of the source tree.
710 723

	
711 724
In order to compile them, use the <tt>make demo</tt> or the
712 725
<tt>make check</tt> commands.
713 726
*/
714 727

	
715 728
}
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
/**
20 20
\mainpage LEMON Documentation
21 21

	
22 22
\section intro Introduction
23 23

	
24
\subsection whatis What is LEMON
25

	
26
LEMON stands for <b>L</b>ibrary for <b>E</b>fficient <b>M</b>odeling
27
and <b>O</b>ptimization in <b>N</b>etworks.
28
It is a C++ template
29
library aimed at combinatorial optimization tasks which
30
often involve in working
31
with graphs.
24
<b>LEMON</b> stands for <i><b>L</b>ibrary for <b>E</b>fficient <b>M</b>odeling
25
and <b>O</b>ptimization in <b>N</b>etworks</i>.
26
It is a C++ template library providing efficient implementation of common
27
data structures and algorithms with focus on combinatorial optimization
28
problems in graphs and networks.
32 29

	
33 30
<b>
34 31
LEMON is an <a class="el" href="http://opensource.org/">open&nbsp;source</a>
35 32
project.
36 33
You are free to use it in your commercial or
37 34
non-commercial applications under very permissive
38 35
\ref license "license terms".
39 36
</b>
40 37

	
41
\subsection howtoread How to read the documentation
38
The project is maintained by the 
39
<a href="http://www.cs.elte.hu/egres/">Egerv&aacute;ry Research Group on
40
Combinatorial Optimization</a> \ref egres
41
at the Operations Research Department of the
42
<a href="http://www.elte.hu/">E&ouml;tv&ouml;s Lor&aacute;nd University,
43
Budapest</a>, Hungary.
44
LEMON is also a member of the <a href="http://www.coin-or.org/">COIN-OR</a>
45
initiative \ref coinor.
46

	
47
\section howtoread How to Read the Documentation
42 48

	
43 49
If you would like to get to know the library, see
44 50
<a class="el" href="http://lemon.cs.elte.hu/pub/tutorial/">LEMON Tutorial</a>.
45 51

	
46 52
If you know what you are looking for, then try to find it under the
47 53
<a class="el" href="modules.html">Modules</a> section.
48 54

	
49 55
If you are a user of the old (0.x) series of LEMON, please check out the
50 56
\ref migration "Migration Guide" for the backward incompatibilities.
51 57
*/
Ignore white space 24576 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
namespace lemon {
20 20

	
21 21
/**
22 22
\page min_cost_flow Minimum Cost Flow Problem
23 23

	
24 24
\section mcf_def Definition (GEQ form)
25 25

	
26 26
The \e minimum \e cost \e flow \e problem is to find a feasible flow of
27 27
minimum total cost from a set of supply nodes to a set of demand nodes
28 28
in a network with capacity constraints (lower and upper bounds)
29
and arc costs.
29
and arc costs \ref amo93networkflows.
30 30

	
31 31
Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$,
32 32
\f$upper: A\rightarrow\mathbf{R}\cup\{+\infty\}\f$ denote the lower and
33 33
upper bounds for the flow values on the arcs, for which
34 34
\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
35 35
\f$cost: A\rightarrow\mathbf{R}\f$ denotes the cost per unit flow
36 36
on the arcs and \f$sup: V\rightarrow\mathbf{R}\f$ denotes the
37 37
signed supply values of the nodes.
38 38
If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
39 39
supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
40 40
\f$-sup(u)\f$ demand.
41 41
A minimum cost flow is an \f$f: A\rightarrow\mathbf{R}\f$ solution
42 42
of the following optimization problem.
43 43

	
44 44
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
45 45
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
46 46
    sup(u) \quad \forall u\in V \f]
47 47
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
48 48

	
49 49
The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
50 50
zero or negative in order to have a feasible solution (since the sum
51 51
of the expressions on the left-hand side of the inequalities is zero).
52 52
It means that the total demand must be greater or equal to the total
53 53
supply and all the supplies have to be carried out from the supply nodes,
54 54
but there could be demands that are not satisfied.
55 55
If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
56 56
constraints have to be satisfied with equality, i.e. all demands
57 57
have to be satisfied and all supplies have to be used.
58 58

	
59 59

	
60 60
\section mcf_algs Algorithms
61 61

	
62 62
LEMON contains several algorithms for solving this problem, for more
63 63
information see \ref min_cost_flow_algs "Minimum Cost Flow Algorithms".
64 64

	
65 65
A feasible solution for this problem can be found using \ref Circulation.
66 66

	
67 67

	
68 68
\section mcf_dual Dual Solution
69 69

	
70 70
The dual solution of the minimum cost flow problem is represented by
71 71
node potentials \f$\pi: V\rightarrow\mathbf{R}\f$.
72 72
An \f$f: A\rightarrow\mathbf{R}\f$ primal feasible solution is optimal
73 73
if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$ node potentials
74 74
the following \e complementary \e slackness optimality conditions hold.
75 75

	
76 76
 - For all \f$uv\in A\f$ arcs:
77 77
   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
78 78
   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
79 79
   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
80 80
 - For all \f$u\in V\f$ nodes:
81 81
   - \f$\pi(u)<=0\f$;
82 82
   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
83 83
     then \f$\pi(u)=0\f$.
84 84
 
85 85
Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
86 86
\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
87 87
\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
88 88

	
89 89
All algorithms provide dual solution (node potentials), as well,
90 90
if an optimal flow is found.
91 91

	
92 92

	
93 93
\section mcf_eq Equality Form
94 94

	
95 95
The above \ref mcf_def "definition" is actually more general than the
96 96
usual formulation of the minimum cost flow problem, in which strict
97 97
equalities are required in the supply/demand contraints.
98 98

	
99 99
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
100 100
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
101 101
    sup(u) \quad \forall u\in V \f]
102 102
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
103 103

	
104 104
However if the sum of the supply values is zero, then these two problems
105 105
are equivalent.
106 106
The \ref min_cost_flow_algs "algorithms" in LEMON support the general
107 107
form, so if you need the equality form, you have to ensure this additional
108 108
contraint manually.
109 109

	
110 110

	
111 111
\section mcf_leq Opposite Inequalites (LEQ Form)
112 112

	
113 113
Another possible definition of the minimum cost flow problem is
114 114
when there are <em>"less or equal"</em> (LEQ) supply/demand constraints,
115 115
instead of the <em>"greater or equal"</em> (GEQ) constraints.
116 116

	
117 117
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
118 118
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
119 119
    sup(u) \quad \forall u\in V \f]
120 120
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
121 121

	
122 122
It means that the total demand must be less or equal to the 
123 123
total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
124 124
positive) and all the demands have to be satisfied, but there
125 125
could be supplies that are not carried out from the supply
126 126
nodes.
127 127
The equality form is also a special case of this form, of course.
128 128

	
129 129
You could easily transform this case to the \ref mcf_def "GEQ form"
130 130
of the problem by reversing the direction of the arcs and taking the
131 131
negative of the supply values (e.g. using \ref ReverseDigraph and
132 132
\ref NegMap adaptors).
133 133
However \ref NetworkSimplex algorithm also supports this form directly
134 134
for the sake of convenience.
135 135

	
136 136
Note that the optimality conditions for this supply constraint type are
137 137
slightly differ from the conditions that are discussed for the GEQ form,
138 138
namely the potentials have to be non-negative instead of non-positive.
139 139
An \f$f: A\rightarrow\mathbf{R}\f$ feasible solution of this problem
140 140
is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$
141 141
node potentials the following conditions hold.
142 142

	
143 143
 - For all \f$uv\in A\f$ arcs:
144 144
   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
145 145
   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
146 146
   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
147 147
 - For all \f$u\in V\f$ nodes:
148 148
   - \f$\pi(u)>=0\f$;
149 149
   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
150 150
     then \f$\pi(u)=0\f$.
151 151

	
152 152
*/
153 153
}
Ignore white space 6 line context
1 1
%%%%% Defining LEMON %%%%%
2 2

	
3 3
@misc{lemon,
4 4
  key =          {LEMON},
5 5
  title =        {{LEMON} -- {L}ibrary for {E}fficient {M}odeling and
6 6
                  {O}ptimization in {N}etworks},
7 7
  howpublished = {\url{http://lemon.cs.elte.hu/}},
8 8
  year =         2009
9 9
}
10 10

	
11 11
@misc{egres,
12 12
  key =          {EGRES},
13 13
  title =        {{EGRES} -- {E}gerv{\'a}ry {R}esearch {G}roup on
14 14
                  {C}ombinatorial {O}ptimization},
15 15
  url =          {http://www.cs.elte.hu/egres/}
16 16
}
17 17

	
18 18
@misc{coinor,
19 19
  key =          {COIN-OR},
20 20
  title =        {{COIN-OR} -- {C}omputational {I}nfrastructure for
21 21
                  {O}perations {R}esearch},
22 22
  url =          {http://www.coin-or.org/}
23 23
}
24 24

	
25 25

	
26 26
%%%%% Other libraries %%%%%%
27 27

	
28 28
@misc{boost,
29 29
  key =          {Boost},
30 30
  title =        {{B}oost {C++} {L}ibraries},
31 31
  url =          {http://www.boost.org/}
32 32
}
33 33

	
34 34
@book{bglbook,
35 35
  author =       {Jeremy G. Siek and Lee-Quan Lee and Andrew
36 36
                  Lumsdaine},
37 37
  title =        {The Boost Graph Library: User Guide and Reference
38 38
                  Manual},
39 39
  publisher =    {Addison-Wesley},
40 40
  year =         2002
41 41
}
42 42

	
43 43
@misc{leda,
44 44
  key =          {LEDA},
45 45
  title =        {{LEDA} -- {L}ibrary of {E}fficient {D}ata {T}ypes and
46 46
                  {A}lgorithms},
47 47
  url =          {http://www.algorithmic-solutions.com/}
48 48
}
49 49

	
50 50
@book{ledabook,
51 51
  author =       {Kurt Mehlhorn and Stefan N{\"a}her},
52 52
  title =        {{LEDA}: {A} platform for combinatorial and geometric
53 53
                  computing},
54 54
  isbn =         {0-521-56329-1},
55 55
  publisher =    {Cambridge University Press},
56 56
  address =      {New York, NY, USA},
57 57
  year =         1999
58 58
}
59 59

	
60 60

	
61 61
%%%%% Tools that LEMON depends on %%%%%
62 62

	
63 63
@misc{cmake,
64 64
  key =          {CMake},
65 65
  title =        {{CMake} -- {C}ross {P}latform {M}ake},
66 66
  url =          {http://www.cmake.org/}
67 67
}
68 68

	
69 69
@misc{doxygen,
70 70
  key =          {Doxygen},
71 71
  title =        {{Doxygen} -- {S}ource code documentation generator
72 72
                  tool},
73 73
  url =          {http://www.doxygen.org/}
74 74
}
75 75

	
76 76

	
77 77
%%%%% LP/MIP libraries %%%%%
78 78

	
79 79
@misc{glpk,
80 80
  key =          {GLPK},
81 81
  title =        {{GLPK} -- {GNU} {L}inear {P}rogramming {K}it},
82 82
  url =          {http://www.gnu.org/software/glpk/}
83 83
}
84 84

	
85 85
@misc{clp,
86 86
  key =          {Clp},
87 87
  title =        {{Clp} -- {Coin-Or} {L}inear {P}rogramming},
88 88
  url =          {http://projects.coin-or.org/Clp/}
89 89
}
90 90

	
91 91
@misc{cbc,
92 92
  key =          {Cbc},
93 93
  title =        {{Cbc} -- {Coin-Or} {B}ranch and {C}ut},
94 94
  url =          {http://projects.coin-or.org/Cbc/}
95 95
}
96 96

	
97 97
@misc{cplex,
98 98
  key =          {CPLEX},
99 99
  title =        {{ILOG} {CPLEX}},
100 100
  url =          {http://www.ilog.com/}
101 101
}
102 102

	
103 103
@misc{soplex,
104 104
  key =          {SoPlex},
105 105
  title =        {{SoPlex} -- {T}he {S}equential {O}bject-{O}riented
106 106
                  {S}implex},
107 107
  url =          {http://soplex.zib.de/}
108 108
}
109 109

	
110 110

	
111 111
%%%%% General books %%%%%
112 112

	
113 113
@book{amo93networkflows,
114 114
  author =       {Ravindra K. Ahuja and Thomas L. Magnanti and James
115 115
                  B. Orlin},
116 116
  title =        {Network Flows: Theory, Algorithms, and Applications},
117 117
  publisher =    {Prentice-Hall, Inc.},
118 118
  year =         1993,
119 119
  month =        feb,
120 120
  isbn =         {978-0136175490}
121 121
}
122 122

	
123 123
@book{schrijver03combinatorial,
124 124
  author =       {Alexander Schrijver},
125 125
  title =        {Combinatorial Optimization: Polyhedra and Efficiency},
126 126
  publisher =    {Springer-Verlag},
127 127
  year =         2003,
128 128
  isbn =         {978-3540443896}
129 129
}
130 130

	
131 131
@book{clrs01algorithms,
132 132
  author =       {Thomas H. Cormen and Charles E. Leiserson and Ronald
133 133
                  L. Rivest and Clifford Stein},
134 134
  title =        {Introduction to Algorithms},
135 135
  publisher =    {The MIT Press},
136 136
  year =         2001,
137 137
  edition =      {2nd}
138 138
}
139 139

	
140 140
@book{stroustrup00cpp,
141 141
  author =       {Bjarne Stroustrup},
142 142
  title =        {The C++ Programming Language},
143 143
  edition =      {3rd},
144 144
  publisher =    {Addison-Wesley Professional},
145 145
  isbn =         0201700735,
146 146
  month =        {February},
147 147
  year =         2000
148 148
}
149 149

	
150 150

	
151 151
%%%%% Maximum flow algorithms %%%%%
152 152

	
153
@inproceedings{goldberg86newapproach,
153
@article{edmondskarp72theoretical,
154
  author =       {Jack Edmonds and Richard M. Karp},
155
  title =        {Theoretical improvements in algorithmic efficiency
156
                  for network flow problems},
157
  journal =      {Journal of the ACM},
158
  year =         1972,
159
  volume =       19,
160
  number =       2,
161
  pages =        {248-264}
162
}
163

	
164
@article{goldberg88newapproach,
154 165
  author =       {Andrew V. Goldberg and Robert E. Tarjan},
155 166
  title =        {A new approach to the maximum flow problem},
156
  booktitle =    {STOC '86: Proceedings of the Eighteenth Annual ACM
157
                  Symposium on Theory of Computing},
158
  year =         1986,
159
  publisher =    {ACM Press},
160
  address =      {New York, NY},
161
  pages =        {136-146}
167
  journal =      {Journal of the ACM},
168
  year =         1988,
169
  volume =       35,
170
  number =       4,
171
  pages =        {921-940}
162 172
}
163 173

	
164 174
@article{dinic70algorithm,
165 175
  author =       {E. A. Dinic},
166 176
  title =        {Algorithm for solution of a problem of maximum flow
167 177
                  in a network with power estimation},
168 178
  journal =      {Soviet Math. Doklady},
169 179
  year =         1970,
170 180
  volume =       11,
171 181
  pages =        {1277-1280}
172 182
}
173 183

	
174 184
@article{goldberg08partial,
175 185
  author =       {Andrew V. Goldberg},
176 186
  title =        {The Partial Augment-Relabel Algorithm for the
177 187
                  Maximum Flow Problem},
178 188
  journal =      {16th Annual European Symposium on Algorithms},
179 189
  year =         2008,
180 190
  pages =        {466-477}
181 191
}
182 192

	
183 193
@article{sleator83dynamic,
184 194
  author =       {Daniel D. Sleator and Robert E. Tarjan},
185 195
  title =        {A data structure for dynamic trees},
186 196
  journal =      {Journal of Computer and System Sciences},
187 197
  year =         1983,
188 198
  volume =       26,
189 199
  number =       3,
190 200
  pages =        {362-391}
191 201
}
192 202

	
193 203

	
194 204
%%%%% Minimum mean cycle algorithms %%%%%
195 205

	
196 206
@article{karp78characterization,
197 207
  author =       {Richard M. Karp},
198 208
  title =        {A characterization of the minimum cycle mean in a
199 209
                  digraph},
200 210
  journal =      {Discrete Math.},
201 211
  year =         1978,
202 212
  volume =       23,
203 213
  pages =        {309-311}
204 214
}
205 215

	
206 216
@article{dasdan98minmeancycle,
207 217
  author =       {Ali Dasdan and Rajesh K. Gupta},
208 218
  title =        {Faster Maximum and Minimum Mean Cycle Alogrithms for
209 219
                  System Performance Analysis},
210 220
  journal =      {IEEE Transactions on Computer-Aided Design of
211 221
                  Integrated Circuits and Systems},
212 222
  year =         1998,
213 223
  volume =       17,
214 224
  number =       10,
215 225
  pages =        {889-899}
216 226
}
217 227

	
218 228

	
219 229
%%%%% Minimum cost flow algorithms %%%%%
220 230

	
221 231
@article{klein67primal,
222 232
  author =       {Morton Klein},
223 233
  title =        {A primal method for minimal cost flows with
224 234
                  applications to the assignment and transportation
225 235
                  problems},
226 236
  journal =      {Management Science},
227 237
  year =         1967,
228 238
  volume =       14,
229 239
  pages =        {205-220}
230 240
}
231 241

	
232
@inproceedings{goldberg88cyclecanceling,
242
@article{goldberg89cyclecanceling,
233 243
  author =       {Andrew V. Goldberg and Robert E. Tarjan},
234 244
  title =        {Finding minimum-cost circulations by canceling
235 245
                  negative cycles},
236
  booktitle =    {STOC '88: Proceedings of the Twentieth Annual ACM
237
                  Symposium on Theory of Computing},
238
  year =         1988,
239
  publisher =    {ACM Press},
240
  address =      {New York, NY},
241
  pages =        {388-397}
246
  journal =      {Journal of the ACM},
247
  year =         1989,
248
  volume =       36,
249
  number =       4,
250
  pages =        {873-886}
242 251
}
243 252

	
244
@article{edmondskarp72theoretical,
245
  author =       {Jack Edmonds and Richard M. Karp},
246
  title =        {Theoretical improvements in algorithmic efficiency
247
                  for network flow problems},
248
  journal =      {Journal of the ACM},
249
  year =         1972,
250
  volume =       19,
251
  number =       2,
252
  pages =        {248-264}
253
}
254

	
255
@inproceedings{goldberg87approximation,
256
  author =       {Andrew V. Goldberg and Robert E. Tarjan},
257
  title =        {Solving minimum-cost flow problems by successive
258
                  approximation},
259
  booktitle =    {STOC '87: Proceedings of the Nineteenth Annual ACM
260
                  Symposium on Theory of Computing},
261
  year =         1987,
262
  publisher =    {ACM Press},
263
  address =      {New York, NY},
264
  pages =        {7-18}
265
}
266

	
267
@article{goldberg90finding,
253
@article{goldberg90approximation,
268 254
  author =       {Andrew V. Goldberg and Robert E. Tarjan},
269 255
  title =        {Finding Minimum-Cost Circulations by Successive
270 256
                  Approximation},
271 257
  journal =      {Mathematics of Operations Research},
272 258
  year =         1990,
273 259
  volume =       15,
274 260
  number =       3,
275 261
  pages =        {430-466}
276 262
}
277 263

	
278 264
@article{goldberg97efficient,
279 265
  author =       {Andrew V. Goldberg},
280 266
  title =        {An Efficient Implementation of a Scaling
281 267
                  Minimum-Cost Flow Algorithm},
282 268
  journal =      {Journal of Algorithms},
283 269
  year =         1997,
284 270
  volume =       22,
285 271
  number =       1,
286 272
  pages =        {1-29}
287 273
}
288 274

	
289 275
@article{bunnagel98efficient,
290 276
  author =       {Ursula B{\"u}nnagel and Bernhard Korte and Jens
291 277
                  Vygen},
292 278
  title =        {Efficient implementation of the {G}oldberg-{T}arjan
293 279
                  minimum-cost flow algorithm},
294 280
  journal =      {Optimization Methods and Software},
295 281
  year =         1998,
296 282
  volume =       10,
297 283
  pages =        {157-174}
298 284
}
299 285

	
286
@book{dantzig63linearprog,
287
  author =       {George B. Dantzig},
288
  title =        {Linear Programming and Extensions},
289
  publisher =    {Princeton University Press},
290
  year =         1963
291
}
292

	
300 293
@mastersthesis{kellyoneill91netsimplex,
301 294
  author =       {Damian J. Kelly and Garrett M. O'Neill},
302 295
  title =        {The Minimum Cost Flow Problem and The Network
303 296
                  Simplex Method},
304 297
  school =       {University College},
305 298
  address =      {Dublin, Ireland},
306 299
  year =         1991,
307 300
  month =        sep,
308 301
}
309

	
310
@techreport{lobel96networksimplex,
311
  author =       {Andreas L{\"o}bel},
312
  title =        {Solving large-scale real-world minimum-cost flow
313
                  problems by a network simplex method},
314
  institution =  {Konrad-Zuse-Zentrum fur Informationstechnik Berlin
315
                  ({ZIB})},
316
  address =      {Berlin, Germany},
317
  year =         1996,
318
  number =       {SC 96-7}
319
}
320

	
321
@article{frangioni06computational,
322
  author =       {Antonio Frangioni and Antonio Manca},
323
  title =        {A Computational Study of Cost Reoptimization for
324
                  Min-Cost Flow Problems},
325
  journal =      {INFORMS Journal On Computing},
326
  year =         2006,
327
  volume =       18,
328
  number =       1,
329
  pages =        {61-70}
330
}
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_NETWORK_SIMPLEX_H
20 20
#define LEMON_NETWORK_SIMPLEX_H
21 21

	
22 22
/// \ingroup min_cost_flow_algs
23 23
///
24 24
/// \file
25 25
/// \brief Network Simplex algorithm for finding a minimum cost flow.
26 26

	
27 27
#include <vector>
28 28
#include <limits>
29 29
#include <algorithm>
30 30

	
31 31
#include <lemon/core.h>
32 32
#include <lemon/math.h>
33 33

	
34 34
namespace lemon {
35 35

	
36 36
  /// \addtogroup min_cost_flow_algs
37 37
  /// @{
38 38

	
39 39
  /// \brief Implementation of the primal Network Simplex algorithm
40 40
  /// for finding a \ref min_cost_flow "minimum cost flow".
41 41
  ///
42 42
  /// \ref NetworkSimplex implements the primal Network Simplex algorithm
43
  /// for finding a \ref min_cost_flow "minimum cost flow".
43
  /// for finding a \ref min_cost_flow "minimum cost flow"
44
  /// \ref amo93networkflows, \ref dantzig63linearprog,
45
  /// \ref kellyoneill91netsimplex.
44 46
  /// This algorithm is a specialized version of the linear programming
45 47
  /// simplex method directly for the minimum cost flow problem.
46 48
  /// It is one of the most efficient solution methods.
47 49
  ///
48 50
  /// In general this class is the fastest implementation available
49 51
  /// in LEMON for the minimum cost flow problem.
50 52
  /// Moreover it supports both directions of the supply/demand inequality
51 53
  /// constraints. For more information see \ref SupplyType.
52 54
  ///
53 55
  /// Most of the parameters of the problem (except for the digraph)
54 56
  /// can be given using separate functions, and the algorithm can be
55 57
  /// executed using the \ref run() function. If some parameters are not
56 58
  /// specified, then default values will be used.
57 59
  ///
58 60
  /// \tparam GR The digraph type the algorithm runs on.
59 61
  /// \tparam V The value type used for flow amounts, capacity bounds
60 62
  /// and supply values in the algorithm. By default it is \c int.
61 63
  /// \tparam C The value type used for costs and potentials in the
62 64
  /// algorithm. By default it is the same as \c V.
63 65
  ///
64 66
  /// \warning Both value types must be signed and all input data must
65 67
  /// be integer.
66 68
  ///
67 69
  /// \note %NetworkSimplex provides five different pivot rule
68 70
  /// implementations, from which the most efficient one is used
69 71
  /// by default. For more information see \ref PivotRule.
70 72
  template <typename GR, typename V = int, typename C = V>
71 73
  class NetworkSimplex
72 74
  {
73 75
  public:
74 76

	
75 77
    /// The type of the flow amounts, capacity bounds and supply values
76 78
    typedef V Value;
77 79
    /// The type of the arc costs
78 80
    typedef C Cost;
79 81

	
80 82
  public:
81 83

	
82 84
    /// \brief Problem type constants for the \c run() function.
83 85
    ///
84 86
    /// Enum type containing the problem type constants that can be
85 87
    /// returned by the \ref run() function of the algorithm.
86 88
    enum ProblemType {
87 89
      /// The problem has no feasible solution (flow).
88 90
      INFEASIBLE,
89 91
      /// The problem has optimal solution (i.e. it is feasible and
90 92
      /// bounded), and the algorithm has found optimal flow and node
91 93
      /// potentials (primal and dual solutions).
92 94
      OPTIMAL,
93 95
      /// The objective function of the problem is unbounded, i.e.
94 96
      /// there is a directed cycle having negative total cost and
95 97
      /// infinite upper bound.
96 98
      UNBOUNDED
97 99
    };
98 100
    
99 101
    /// \brief Constants for selecting the type of the supply constraints.
100 102
    ///
101 103
    /// Enum type containing constants for selecting the supply type,
102 104
    /// i.e. the direction of the inequalities in the supply/demand
103 105
    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
104 106
    ///
105 107
    /// The default supply type is \c GEQ, the \c LEQ type can be
106 108
    /// selected using \ref supplyType().
107 109
    /// The equality form is a special case of both supply types.
108 110
    enum SupplyType {
109 111
      /// This option means that there are <em>"greater or equal"</em>
110 112
      /// supply/demand constraints in the definition of the problem.
111 113
      GEQ,
112 114
      /// This option means that there are <em>"less or equal"</em>
113 115
      /// supply/demand constraints in the definition of the problem.
114 116
      LEQ
115 117
    };
116 118
    
117 119
    /// \brief Constants for selecting the pivot rule.
118 120
    ///
119 121
    /// Enum type containing constants for selecting the pivot rule for
120 122
    /// the \ref run() function.
121 123
    ///
122 124
    /// \ref NetworkSimplex provides five different pivot rule
123 125
    /// implementations that significantly affect the running time
124 126
    /// of the algorithm.
125 127
    /// By default \ref BLOCK_SEARCH "Block Search" is used, which
126 128
    /// proved to be the most efficient and the most robust on various
127 129
    /// test inputs according to our benchmark tests.
128 130
    /// However another pivot rule can be selected using the \ref run()
129 131
    /// function with the proper parameter.
130 132
    enum PivotRule {
131 133

	
132 134
      /// The First Eligible pivot rule.
133 135
      /// The next eligible arc is selected in a wraparound fashion
134 136
      /// in every iteration.
135 137
      FIRST_ELIGIBLE,
136 138

	
137 139
      /// The Best Eligible pivot rule.
138 140
      /// The best eligible arc is selected in every iteration.
139 141
      BEST_ELIGIBLE,
140 142

	
141 143
      /// The Block Search pivot rule.
142 144
      /// A specified number of arcs are examined in every iteration
143 145
      /// in a wraparound fashion and the best eligible arc is selected
144 146
      /// from this block.
145 147
      BLOCK_SEARCH,
146 148

	
147 149
      /// The Candidate List pivot rule.
148 150
      /// In a major iteration a candidate list is built from eligible arcs
149 151
      /// in a wraparound fashion and in the following minor iterations
150 152
      /// the best eligible arc is selected from this list.
151 153
      CANDIDATE_LIST,
152 154

	
153 155
      /// The Altering Candidate List pivot rule.
154 156
      /// It is a modified version of the Candidate List method.
155 157
      /// It keeps only the several best eligible arcs from the former
156 158
      /// candidate list and extends this list in every iteration.
157 159
      ALTERING_LIST
158 160
    };
159 161
    
160 162
  private:
161 163

	
162 164
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
163 165

	
164 166
    typedef std::vector<int> IntVector;
165 167
    typedef std::vector<bool> BoolVector;
166 168
    typedef std::vector<Value> ValueVector;
167 169
    typedef std::vector<Cost> CostVector;
168 170

	
169 171
    // State constants for arcs
170 172
    enum ArcStateEnum {
171 173
      STATE_UPPER = -1,
172 174
      STATE_TREE  =  0,
173 175
      STATE_LOWER =  1
174 176
    };
175 177

	
176 178
  private:
177 179

	
178 180
    // Data related to the underlying digraph
179 181
    const GR &_graph;
180 182
    int _node_num;
181 183
    int _arc_num;
182 184
    int _all_arc_num;
183 185
    int _search_arc_num;
184 186

	
185 187
    // Parameters of the problem
186 188
    bool _have_lower;
187 189
    SupplyType _stype;
188 190
    Value _sum_supply;
189 191

	
190 192
    // Data structures for storing the digraph
191 193
    IntNodeMap _node_id;
192 194
    IntArcMap _arc_id;
193 195
    IntVector _source;
194 196
    IntVector _target;
195 197

	
196 198
    // Node and arc data
197 199
    ValueVector _lower;
198 200
    ValueVector _upper;
199 201
    ValueVector _cap;
200 202
    CostVector _cost;
201 203
    ValueVector _supply;
202 204
    ValueVector _flow;
203 205
    CostVector _pi;
204 206

	
205 207
    // Data for storing the spanning tree structure
206 208
    IntVector _parent;
207 209
    IntVector _pred;
208 210
    IntVector _thread;
209 211
    IntVector _rev_thread;
210 212
    IntVector _succ_num;
211 213
    IntVector _last_succ;
212 214
    IntVector _dirty_revs;
213 215
    BoolVector _forward;
214 216
    IntVector _state;
215 217
    int _root;
216 218

	
217 219
    // Temporary data used in the current pivot iteration
218 220
    int in_arc, join, u_in, v_in, u_out, v_out;
219 221
    int first, second, right, last;
220 222
    int stem, par_stem, new_stem;
221 223
    Value delta;
222 224

	
223 225
  public:
224 226
  
225 227
    /// \brief Constant for infinite upper bounds (capacities).
226 228
    ///
227 229
    /// Constant for infinite upper bounds (capacities).
228 230
    /// It is \c std::numeric_limits<Value>::infinity() if available,
229 231
    /// \c std::numeric_limits<Value>::max() otherwise.
230 232
    const Value INF;
231 233

	
232 234
  private:
233 235

	
234 236
    // Implementation of the First Eligible pivot rule
235 237
    class FirstEligiblePivotRule
236 238
    {
237 239
    private:
238 240

	
239 241
      // References to the NetworkSimplex class
240 242
      const IntVector  &_source;
241 243
      const IntVector  &_target;
242 244
      const CostVector &_cost;
243 245
      const IntVector  &_state;
244 246
      const CostVector &_pi;
245 247
      int &_in_arc;
246 248
      int _search_arc_num;
247 249

	
248 250
      // Pivot rule data
249 251
      int _next_arc;
250 252

	
251 253
    public:
252 254

	
253 255
      // Constructor
254 256
      FirstEligiblePivotRule(NetworkSimplex &ns) :
255 257
        _source(ns._source), _target(ns._target),
256 258
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
257 259
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
258 260
        _next_arc(0)
259 261
      {}
260 262

	
261 263
      // Find next entering arc
262 264
      bool findEnteringArc() {
263 265
        Cost c;
264 266
        for (int e = _next_arc; e < _search_arc_num; ++e) {
265 267
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
266 268
          if (c < 0) {
267 269
            _in_arc = e;
268 270
            _next_arc = e + 1;
269 271
            return true;
270 272
          }
271 273
        }
272 274
        for (int e = 0; e < _next_arc; ++e) {
273 275
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
274 276
          if (c < 0) {
275 277
            _in_arc = e;
276 278
            _next_arc = e + 1;
277 279
            return true;
278 280
          }
279 281
        }
280 282
        return false;
281 283
      }
282 284

	
283 285
    }; //class FirstEligiblePivotRule
284 286

	
285 287

	
286 288
    // Implementation of the Best Eligible pivot rule
287 289
    class BestEligiblePivotRule
288 290
    {
289 291
    private:
290 292

	
291 293
      // References to the NetworkSimplex class
292 294
      const IntVector  &_source;
293 295
      const IntVector  &_target;
294 296
      const CostVector &_cost;
295 297
      const IntVector  &_state;
296 298
      const CostVector &_pi;
297 299
      int &_in_arc;
298 300
      int _search_arc_num;
299 301

	
300 302
    public:
301 303

	
302 304
      // Constructor
303 305
      BestEligiblePivotRule(NetworkSimplex &ns) :
304 306
        _source(ns._source), _target(ns._target),
305 307
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
306 308
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
307 309
      {}
308 310

	
309 311
      // Find next entering arc
310 312
      bool findEnteringArc() {
311 313
        Cost c, min = 0;
312 314
        for (int e = 0; e < _search_arc_num; ++e) {
313 315
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
314 316
          if (c < min) {
315 317
            min = c;
316 318
            _in_arc = e;
317 319
          }
318 320
        }
319 321
        return min < 0;
320 322
      }
321 323

	
322 324
    }; //class BestEligiblePivotRule
323 325

	
324 326

	
325 327
    // Implementation of the Block Search pivot rule
326 328
    class BlockSearchPivotRule
327 329
    {
328 330
    private:
329 331

	
330 332
      // References to the NetworkSimplex class
331 333
      const IntVector  &_source;
332 334
      const IntVector  &_target;
333 335
      const CostVector &_cost;
334 336
      const IntVector  &_state;
335 337
      const CostVector &_pi;
336 338
      int &_in_arc;
337 339
      int _search_arc_num;
338 340

	
339 341
      // Pivot rule data
340 342
      int _block_size;
341 343
      int _next_arc;
342 344

	
343 345
    public:
344 346

	
345 347
      // Constructor
346 348
      BlockSearchPivotRule(NetworkSimplex &ns) :
347 349
        _source(ns._source), _target(ns._target),
348 350
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
349 351
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
350 352
        _next_arc(0)
351 353
      {
352 354
        // The main parameters of the pivot rule
353 355
        const double BLOCK_SIZE_FACTOR = 0.5;
354 356
        const int MIN_BLOCK_SIZE = 10;
355 357

	
356 358
        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
357 359
                                    std::sqrt(double(_search_arc_num))),
358 360
                                MIN_BLOCK_SIZE );
359 361
      }
360 362

	
361 363
      // Find next entering arc
362 364
      bool findEnteringArc() {
363 365
        Cost c, min = 0;
364 366
        int cnt = _block_size;
365 367
        int e;
366 368
        for (e = _next_arc; e < _search_arc_num; ++e) {
367 369
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
368 370
          if (c < min) {
369 371
            min = c;
370 372
            _in_arc = e;
371 373
          }
372 374
          if (--cnt == 0) {
373 375
            if (min < 0) goto search_end;
374 376
            cnt = _block_size;
375 377
          }
376 378
        }
377 379
        for (e = 0; e < _next_arc; ++e) {
378 380
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
379 381
          if (c < min) {
380 382
            min = c;
381 383
            _in_arc = e;
382 384
          }
383 385
          if (--cnt == 0) {
384 386
            if (min < 0) goto search_end;
385 387
            cnt = _block_size;
386 388
          }
387 389
        }
388 390
        if (min >= 0) return false;
389 391

	
390 392
      search_end:
391 393
        _next_arc = e;
392 394
        return true;
393 395
      }
394 396

	
395 397
    }; //class BlockSearchPivotRule
396 398

	
397 399

	
398 400
    // Implementation of the Candidate List pivot rule
399 401
    class CandidateListPivotRule
400 402
    {
401 403
    private:
402 404

	
403 405
      // References to the NetworkSimplex class
404 406
      const IntVector  &_source;
405 407
      const IntVector  &_target;
406 408
      const CostVector &_cost;
407 409
      const IntVector  &_state;
408 410
      const CostVector &_pi;
409 411
      int &_in_arc;
410 412
      int _search_arc_num;
411 413

	
412 414
      // Pivot rule data
413 415
      IntVector _candidates;
414 416
      int _list_length, _minor_limit;
415 417
      int _curr_length, _minor_count;
416 418
      int _next_arc;
417 419

	
418 420
    public:
419 421

	
420 422
      /// Constructor
421 423
      CandidateListPivotRule(NetworkSimplex &ns) :
422 424
        _source(ns._source), _target(ns._target),
423 425
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
424 426
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
425 427
        _next_arc(0)
426 428
      {
427 429
        // The main parameters of the pivot rule
428 430
        const double LIST_LENGTH_FACTOR = 0.25;
429 431
        const int MIN_LIST_LENGTH = 10;
430 432
        const double MINOR_LIMIT_FACTOR = 0.1;
431 433
        const int MIN_MINOR_LIMIT = 3;
432 434

	
433 435
        _list_length = std::max( int(LIST_LENGTH_FACTOR *
434 436
                                     std::sqrt(double(_search_arc_num))),
435 437
                                 MIN_LIST_LENGTH );
436 438
        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
437 439
                                 MIN_MINOR_LIMIT );
438 440
        _curr_length = _minor_count = 0;
439 441
        _candidates.resize(_list_length);
440 442
      }
441 443

	
442 444
      /// Find next entering arc
443 445
      bool findEnteringArc() {
444 446
        Cost min, c;
445 447
        int e;
446 448
        if (_curr_length > 0 && _minor_count < _minor_limit) {
447 449
          // Minor iteration: select the best eligible arc from the
448 450
          // current candidate list
449 451
          ++_minor_count;
450 452
          min = 0;
451 453
          for (int i = 0; i < _curr_length; ++i) {
452 454
            e = _candidates[i];
453 455
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
454 456
            if (c < min) {
455 457
              min = c;
456 458
              _in_arc = e;
457 459
            }
458 460
            else if (c >= 0) {
459 461
              _candidates[i--] = _candidates[--_curr_length];
460 462
            }
461 463
          }
462 464
          if (min < 0) return true;
463 465
        }
464 466

	
465 467
        // Major iteration: build a new candidate list
466 468
        min = 0;
467 469
        _curr_length = 0;
468 470
        for (e = _next_arc; e < _search_arc_num; ++e) {
469 471
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
470 472
          if (c < 0) {
471 473
            _candidates[_curr_length++] = e;
472 474
            if (c < min) {
473 475
              min = c;
474 476
              _in_arc = e;
475 477
            }
476 478
            if (_curr_length == _list_length) goto search_end;
477 479
          }
478 480
        }
479 481
        for (e = 0; e < _next_arc; ++e) {
480 482
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
481 483
          if (c < 0) {
482 484
            _candidates[_curr_length++] = e;
483 485
            if (c < min) {
484 486
              min = c;
485 487
              _in_arc = e;
486 488
            }
487 489
            if (_curr_length == _list_length) goto search_end;
488 490
          }
489 491
        }
490 492
        if (_curr_length == 0) return false;
491 493
      
492 494
      search_end:        
493 495
        _minor_count = 1;
494 496
        _next_arc = e;
495 497
        return true;
496 498
      }
497 499

	
498 500
    }; //class CandidateListPivotRule
499 501

	
500 502

	
501 503
    // Implementation of the Altering Candidate List pivot rule
502 504
    class AlteringListPivotRule
503 505
    {
504 506
    private:
505 507

	
506 508
      // References to the NetworkSimplex class
507 509
      const IntVector  &_source;
508 510
      const IntVector  &_target;
509 511
      const CostVector &_cost;
510 512
      const IntVector  &_state;
511 513
      const CostVector &_pi;
512 514
      int &_in_arc;
513 515
      int _search_arc_num;
514 516

	
515 517
      // Pivot rule data
516 518
      int _block_size, _head_length, _curr_length;
517 519
      int _next_arc;
518 520
      IntVector _candidates;
519 521
      CostVector _cand_cost;
520 522

	
521 523
      // Functor class to compare arcs during sort of the candidate list
522 524
      class SortFunc
523 525
      {
524 526
      private:
525 527
        const CostVector &_map;
526 528
      public:
527 529
        SortFunc(const CostVector &map) : _map(map) {}
528 530
        bool operator()(int left, int right) {
529 531
          return _map[left] > _map[right];
530 532
        }
531 533
      };
532 534

	
533 535
      SortFunc _sort_func;
534 536

	
535 537
    public:
536 538

	
537 539
      // Constructor
538 540
      AlteringListPivotRule(NetworkSimplex &ns) :
539 541
        _source(ns._source), _target(ns._target),
540 542
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
541 543
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
542 544
        _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
543 545
      {
544 546
        // The main parameters of the pivot rule
545 547
        const double BLOCK_SIZE_FACTOR = 1.0;
546 548
        const int MIN_BLOCK_SIZE = 10;
547 549
        const double HEAD_LENGTH_FACTOR = 0.1;
548 550
        const int MIN_HEAD_LENGTH = 3;
549 551

	
550 552
        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
551 553
                                    std::sqrt(double(_search_arc_num))),
552 554
                                MIN_BLOCK_SIZE );
553 555
        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
554 556
                                 MIN_HEAD_LENGTH );
555 557
        _candidates.resize(_head_length + _block_size);
556 558
        _curr_length = 0;
557 559
      }
558 560

	
559 561
      // Find next entering arc
560 562
      bool findEnteringArc() {
561 563
        // Check the current candidate list
562 564
        int e;
563 565
        for (int i = 0; i < _curr_length; ++i) {
564 566
          e = _candidates[i];
565 567
          _cand_cost[e] = _state[e] *
566 568
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
567 569
          if (_cand_cost[e] >= 0) {
568 570
            _candidates[i--] = _candidates[--_curr_length];
569 571
          }
570 572
        }
571 573

	
572 574
        // Extend the list
573 575
        int cnt = _block_size;
574 576
        int limit = _head_length;
575 577

	
576 578
        for (e = _next_arc; e < _search_arc_num; ++e) {
577 579
          _cand_cost[e] = _state[e] *
578 580
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
579 581
          if (_cand_cost[e] < 0) {
580 582
            _candidates[_curr_length++] = e;
581 583
          }
582 584
          if (--cnt == 0) {
583 585
            if (_curr_length > limit) goto search_end;
584 586
            limit = 0;
585 587
            cnt = _block_size;
586 588
          }
587 589
        }
588 590
        for (e = 0; e < _next_arc; ++e) {
589 591
          _cand_cost[e] = _state[e] *
590 592
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
591 593
          if (_cand_cost[e] < 0) {
592 594
            _candidates[_curr_length++] = e;
593 595
          }
594 596
          if (--cnt == 0) {
595 597
            if (_curr_length > limit) goto search_end;
596 598
            limit = 0;
597 599
            cnt = _block_size;
598 600
          }
599 601
        }
600 602
        if (_curr_length == 0) return false;
601 603
        
602 604
      search_end:
603 605

	
604 606
        // Make heap of the candidate list (approximating a partial sort)
605 607
        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
606 608
                   _sort_func );
607 609

	
608 610
        // Pop the first element of the heap
609 611
        _in_arc = _candidates[0];
610 612
        _next_arc = e;
611 613
        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
612 614
                  _sort_func );
613 615
        _curr_length = std::min(_head_length, _curr_length - 1);
614 616
        return true;
615 617
      }
616 618

	
617 619
    }; //class AlteringListPivotRule
618 620

	
619 621
  public:
620 622

	
621 623
    /// \brief Constructor.
622 624
    ///
623 625
    /// The constructor of the class.
624 626
    ///
625 627
    /// \param graph The digraph the algorithm runs on.
626 628
    /// \param arc_mixing Indicate if the arcs have to be stored in a
627 629
    /// mixed order in the internal data structure. 
628 630
    /// In special cases, it could lead to better overall performance,
629 631
    /// but it is usually slower. Therefore it is disabled by default.
630 632
    NetworkSimplex(const GR& graph, bool arc_mixing = false) :
631 633
      _graph(graph), _node_id(graph), _arc_id(graph),
632 634
      INF(std::numeric_limits<Value>::has_infinity ?
633 635
          std::numeric_limits<Value>::infinity() :
634 636
          std::numeric_limits<Value>::max())
635 637
    {
636 638
      // Check the value types
637 639
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
638 640
        "The flow type of NetworkSimplex must be signed");
639 641
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
640 642
        "The cost type of NetworkSimplex must be signed");
641 643
        
642 644
      // Resize vectors
643 645
      _node_num = countNodes(_graph);
644 646
      _arc_num = countArcs(_graph);
645 647
      int all_node_num = _node_num + 1;
646 648
      int max_arc_num = _arc_num + 2 * _node_num;
647 649

	
648 650
      _source.resize(max_arc_num);
649 651
      _target.resize(max_arc_num);
650 652

	
651 653
      _lower.resize(_arc_num);
652 654
      _upper.resize(_arc_num);
653 655
      _cap.resize(max_arc_num);
654 656
      _cost.resize(max_arc_num);
655 657
      _supply.resize(all_node_num);
656 658
      _flow.resize(max_arc_num);
657 659
      _pi.resize(all_node_num);
658 660

	
659 661
      _parent.resize(all_node_num);
660 662
      _pred.resize(all_node_num);
661 663
      _forward.resize(all_node_num);
662 664
      _thread.resize(all_node_num);
663 665
      _rev_thread.resize(all_node_num);
664 666
      _succ_num.resize(all_node_num);
665 667
      _last_succ.resize(all_node_num);
666 668
      _state.resize(max_arc_num);
667 669

	
668 670
      // Copy the graph
669 671
      int i = 0;
670 672
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
671 673
        _node_id[n] = i;
672 674
      }
673 675
      if (arc_mixing) {
674 676
        // Store the arcs in a mixed order
675 677
        int k = std::max(int(std::sqrt(double(_arc_num))), 10);
676 678
        int i = 0, j = 0;
677 679
        for (ArcIt a(_graph); a != INVALID; ++a) {
678 680
          _arc_id[a] = i;
679 681
          _source[i] = _node_id[_graph.source(a)];
680 682
          _target[i] = _node_id[_graph.target(a)];
681 683
          if ((i += k) >= _arc_num) i = ++j;
682 684
        }
683 685
      } else {
684 686
        // Store the arcs in the original order
685 687
        int i = 0;
686 688
        for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
687 689
          _arc_id[a] = i;
688 690
          _source[i] = _node_id[_graph.source(a)];
689 691
          _target[i] = _node_id[_graph.target(a)];
690 692
        }
691 693
      }
692 694
      
693 695
      // Reset parameters
694 696
      reset();
695 697
    }
696 698

	
697 699
    /// \name Parameters
698 700
    /// The parameters of the algorithm can be specified using these
699 701
    /// functions.
700 702

	
701 703
    /// @{
702 704

	
703 705
    /// \brief Set the lower bounds on the arcs.
704 706
    ///
705 707
    /// This function sets the lower bounds on the arcs.
706 708
    /// If it is not used before calling \ref run(), the lower bounds
707 709
    /// will be set to zero on all arcs.
708 710
    ///
709 711
    /// \param map An arc map storing the lower bounds.
710 712
    /// Its \c Value type must be convertible to the \c Value type
711 713
    /// of the algorithm.
712 714
    ///
713 715
    /// \return <tt>(*this)</tt>
714 716
    template <typename LowerMap>
715 717
    NetworkSimplex& lowerMap(const LowerMap& map) {
716 718
      _have_lower = true;
717 719
      for (ArcIt a(_graph); a != INVALID; ++a) {
718 720
        _lower[_arc_id[a]] = map[a];
719 721
      }
720 722
      return *this;
721 723
    }
722 724

	
723 725
    /// \brief Set the upper bounds (capacities) on the arcs.
724 726
    ///
725 727
    /// This function sets the upper bounds (capacities) on the arcs.
726 728
    /// If it is not used before calling \ref run(), the upper bounds
727 729
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
728 730
    /// unbounded from above on each arc).
729 731
    ///
730 732
    /// \param map An arc map storing the upper bounds.
731 733
    /// Its \c Value type must be convertible to the \c Value type
732 734
    /// of the algorithm.
733 735
    ///
734 736
    /// \return <tt>(*this)</tt>
735 737
    template<typename UpperMap>
736 738
    NetworkSimplex& upperMap(const UpperMap& map) {
737 739
      for (ArcIt a(_graph); a != INVALID; ++a) {
738 740
        _upper[_arc_id[a]] = map[a];
739 741
      }
740 742
      return *this;
741 743
    }
742 744

	
743 745
    /// \brief Set the costs of the arcs.
744 746
    ///
745 747
    /// This function sets the costs of the arcs.
746 748
    /// If it is not used before calling \ref run(), the costs
747 749
    /// will be set to \c 1 on all arcs.
748 750
    ///
749 751
    /// \param map An arc map storing the costs.
750 752
    /// Its \c Value type must be convertible to the \c Cost type
751 753
    /// of the algorithm.
752 754
    ///
753 755
    /// \return <tt>(*this)</tt>
754 756
    template<typename CostMap>
755 757
    NetworkSimplex& costMap(const CostMap& map) {
756 758
      for (ArcIt a(_graph); a != INVALID; ++a) {
757 759
        _cost[_arc_id[a]] = map[a];
758 760
      }
759 761
      return *this;
760 762
    }
761 763

	
762 764
    /// \brief Set the supply values of the nodes.
763 765
    ///
764 766
    /// This function sets the supply values of the nodes.
765 767
    /// If neither this function nor \ref stSupply() is used before
766 768
    /// calling \ref run(), the supply of each node will be set to zero.
767 769
    ///
768 770
    /// \param map A node map storing the supply values.
769 771
    /// Its \c Value type must be convertible to the \c Value type
770 772
    /// of the algorithm.
771 773
    ///
772 774
    /// \return <tt>(*this)</tt>
773 775
    template<typename SupplyMap>
774 776
    NetworkSimplex& supplyMap(const SupplyMap& map) {
775 777
      for (NodeIt n(_graph); n != INVALID; ++n) {
776 778
        _supply[_node_id[n]] = map[n];
777 779
      }
778 780
      return *this;
779 781
    }
780 782

	
781 783
    /// \brief Set single source and target nodes and a supply value.
782 784
    ///
783 785
    /// This function sets a single source node and a single target node
784 786
    /// and the required flow value.
785 787
    /// If neither this function nor \ref supplyMap() is used before
786 788
    /// calling \ref run(), the supply of each node will be set to zero.
787 789
    ///
788 790
    /// Using this function has the same effect as using \ref supplyMap()
789 791
    /// with such a map in which \c k is assigned to \c s, \c -k is
790 792
    /// assigned to \c t and all other nodes have zero supply value.
791 793
    ///
792 794
    /// \param s The source node.
793 795
    /// \param t The target node.
794 796
    /// \param k The required amount of flow from node \c s to node \c t
795 797
    /// (i.e. the supply of \c s and the demand of \c t).
796 798
    ///
797 799
    /// \return <tt>(*this)</tt>
798 800
    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
799 801
      for (int i = 0; i != _node_num; ++i) {
800 802
        _supply[i] = 0;
801 803
      }
802 804
      _supply[_node_id[s]] =  k;
803 805
      _supply[_node_id[t]] = -k;
804 806
      return *this;
805 807
    }
806 808
    
807 809
    /// \brief Set the type of the supply constraints.
808 810
    ///
809 811
    /// This function sets the type of the supply/demand constraints.
810 812
    /// If it is not used before calling \ref run(), the \ref GEQ supply
811 813
    /// type will be used.
812 814
    ///
813 815
    /// For more information see \ref SupplyType.
814 816
    ///
815 817
    /// \return <tt>(*this)</tt>
816 818
    NetworkSimplex& supplyType(SupplyType supply_type) {
817 819
      _stype = supply_type;
818 820
      return *this;
819 821
    }
820 822

	
821 823
    /// @}
822 824

	
823 825
    /// \name Execution Control
824 826
    /// The algorithm can be executed using \ref run().
825 827

	
826 828
    /// @{
827 829

	
828 830
    /// \brief Run the algorithm.
829 831
    ///
830 832
    /// This function runs the algorithm.
831 833
    /// The paramters can be specified using functions \ref lowerMap(),
832 834
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 
833 835
    /// \ref supplyType().
834 836
    /// For example,
835 837
    /// \code
836 838
    ///   NetworkSimplex<ListDigraph> ns(graph);
837 839
    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
838 840
    ///     .supplyMap(sup).run();
839 841
    /// \endcode
840 842
    ///
841 843
    /// This function can be called more than once. All the parameters
842 844
    /// that have been given are kept for the next call, unless
843 845
    /// \ref reset() is called, thus only the modified parameters
844 846
    /// have to be set again. See \ref reset() for examples.
845 847
    /// However the underlying digraph must not be modified after this
846 848
    /// class have been constructed, since it copies and extends the graph.
847 849
    ///
848 850
    /// \param pivot_rule The pivot rule that will be used during the
849 851
    /// algorithm. For more information see \ref PivotRule.
850 852
    ///
851 853
    /// \return \c INFEASIBLE if no feasible flow exists,
852 854
    /// \n \c OPTIMAL if the problem has optimal solution
853 855
    /// (i.e. it is feasible and bounded), and the algorithm has found
854 856
    /// optimal flow and node potentials (primal and dual solutions),
855 857
    /// \n \c UNBOUNDED if the objective function of the problem is
856 858
    /// unbounded, i.e. there is a directed cycle having negative total
857 859
    /// cost and infinite upper bound.
858 860
    ///
859 861
    /// \see ProblemType, PivotRule
860 862
    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
861 863
      if (!init()) return INFEASIBLE;
862 864
      return start(pivot_rule);
863 865
    }
864 866

	
865 867
    /// \brief Reset all the parameters that have been given before.
866 868
    ///
867 869
    /// This function resets all the paramaters that have been given
868 870
    /// before using functions \ref lowerMap(), \ref upperMap(),
869 871
    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
870 872
    ///
871 873
    /// It is useful for multiple run() calls. If this function is not
872 874
    /// used, all the parameters given before are kept for the next
873 875
    /// \ref run() call.
874 876
    /// However the underlying digraph must not be modified after this
875 877
    /// class have been constructed, since it copies and extends the graph.
876 878
    ///
877 879
    /// For example,
878 880
    /// \code
879 881
    ///   NetworkSimplex<ListDigraph> ns(graph);
880 882
    ///
881 883
    ///   // First run
882 884
    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
883 885
    ///     .supplyMap(sup).run();
884 886
    ///
885 887
    ///   // Run again with modified cost map (reset() is not called,
886 888
    ///   // so only the cost map have to be set again)
887 889
    ///   cost[e] += 100;
888 890
    ///   ns.costMap(cost).run();
889 891
    ///
890 892
    ///   // Run again from scratch using reset()
891 893
    ///   // (the lower bounds will be set to zero on all arcs)
892 894
    ///   ns.reset();
893 895
    ///   ns.upperMap(capacity).costMap(cost)
894 896
    ///     .supplyMap(sup).run();
895 897
    /// \endcode
896 898
    ///
897 899
    /// \return <tt>(*this)</tt>
898 900
    NetworkSimplex& reset() {
899 901
      for (int i = 0; i != _node_num; ++i) {
900 902
        _supply[i] = 0;
901 903
      }
902 904
      for (int i = 0; i != _arc_num; ++i) {
903 905
        _lower[i] = 0;
904 906
        _upper[i] = INF;
905 907
        _cost[i] = 1;
906 908
      }
907 909
      _have_lower = false;
908 910
      _stype = GEQ;
909 911
      return *this;
910 912
    }
911 913

	
912 914
    /// @}
913 915

	
914 916
    /// \name Query Functions
915 917
    /// The results of the algorithm can be obtained using these
916 918
    /// functions.\n
917 919
    /// The \ref run() function must be called before using them.
918 920

	
919 921
    /// @{
920 922

	
921 923
    /// \brief Return the total cost of the found flow.
922 924
    ///
923 925
    /// This function returns the total cost of the found flow.
924 926
    /// Its complexity is O(e).
925 927
    ///
926 928
    /// \note The return type of the function can be specified as a
927 929
    /// template parameter. For example,
928 930
    /// \code
929 931
    ///   ns.totalCost<double>();
930 932
    /// \endcode
931 933
    /// It is useful if the total cost cannot be stored in the \c Cost
932 934
    /// type of the algorithm, which is the default return type of the
933 935
    /// function.
934 936
    ///
935 937
    /// \pre \ref run() must be called before using this function.
936 938
    template <typename Number>
937 939
    Number totalCost() const {
938 940
      Number c = 0;
939 941
      for (ArcIt a(_graph); a != INVALID; ++a) {
940 942
        int i = _arc_id[a];
941 943
        c += Number(_flow[i]) * Number(_cost[i]);
942 944
      }
943 945
      return c;
944 946
    }
945 947

	
946 948
#ifndef DOXYGEN
947 949
    Cost totalCost() const {
948 950
      return totalCost<Cost>();
949 951
    }
950 952
#endif
951 953

	
952 954
    /// \brief Return the flow on the given arc.
953 955
    ///
954 956
    /// This function returns the flow on the given arc.
955 957
    ///
956 958
    /// \pre \ref run() must be called before using this function.
957 959
    Value flow(const Arc& a) const {
958 960
      return _flow[_arc_id[a]];
959 961
    }
960 962

	
961 963
    /// \brief Return the flow map (the primal solution).
962 964
    ///
963 965
    /// This function copies the flow value on each arc into the given
964 966
    /// map. The \c Value type of the algorithm must be convertible to
965 967
    /// the \c Value type of the map.
966 968
    ///
967 969
    /// \pre \ref run() must be called before using this function.
968 970
    template <typename FlowMap>
969 971
    void flowMap(FlowMap &map) const {
970 972
      for (ArcIt a(_graph); a != INVALID; ++a) {
971 973
        map.set(a, _flow[_arc_id[a]]);
972 974
      }
973 975
    }
974 976

	
975 977
    /// \brief Return the potential (dual value) of the given node.
976 978
    ///
977 979
    /// This function returns the potential (dual value) of the
978 980
    /// given node.
979 981
    ///
980 982
    /// \pre \ref run() must be called before using this function.
981 983
    Cost potential(const Node& n) const {
982 984
      return _pi[_node_id[n]];
983 985
    }
984 986

	
985 987
    /// \brief Return the potential map (the dual solution).
986 988
    ///
987 989
    /// This function copies the potential (dual value) of each node
988 990
    /// into the given map.
989 991
    /// The \c Cost type of the algorithm must be convertible to the
990 992
    /// \c Value type of the map.
991 993
    ///
992 994
    /// \pre \ref run() must be called before using this function.
993 995
    template <typename PotentialMap>
994 996
    void potentialMap(PotentialMap &map) const {
995 997
      for (NodeIt n(_graph); n != INVALID; ++n) {
996 998
        map.set(n, _pi[_node_id[n]]);
997 999
      }
998 1000
    }
999 1001

	
1000 1002
    /// @}
1001 1003

	
1002 1004
  private:
1003 1005

	
1004 1006
    // Initialize internal data structures
1005 1007
    bool init() {
1006 1008
      if (_node_num == 0) return false;
1007 1009

	
1008 1010
      // Check the sum of supply values
1009 1011
      _sum_supply = 0;
1010 1012
      for (int i = 0; i != _node_num; ++i) {
1011 1013
        _sum_supply += _supply[i];
1012 1014
      }
1013 1015
      if ( !((_stype == GEQ && _sum_supply <= 0) ||
1014 1016
             (_stype == LEQ && _sum_supply >= 0)) ) return false;
1015 1017

	
1016 1018
      // Remove non-zero lower bounds
1017 1019
      if (_have_lower) {
1018 1020
        for (int i = 0; i != _arc_num; ++i) {
1019 1021
          Value c = _lower[i];
1020 1022
          if (c >= 0) {
1021 1023
            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
1022 1024
          } else {
1023 1025
            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
1024 1026
          }
1025 1027
          _supply[_source[i]] -= c;
1026 1028
          _supply[_target[i]] += c;
1027 1029
        }
1028 1030
      } else {
1029 1031
        for (int i = 0; i != _arc_num; ++i) {
1030 1032
          _cap[i] = _upper[i];
1031 1033
        }
1032 1034
      }
1033 1035

	
1034 1036
      // Initialize artifical cost
1035 1037
      Cost ART_COST;
1036 1038
      if (std::numeric_limits<Cost>::is_exact) {
1037 1039
        ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
1038 1040
      } else {
1039 1041
        ART_COST = std::numeric_limits<Cost>::min();
1040 1042
        for (int i = 0; i != _arc_num; ++i) {
1041 1043
          if (_cost[i] > ART_COST) ART_COST = _cost[i];
1042 1044
        }
1043 1045
        ART_COST = (ART_COST + 1) * _node_num;
1044 1046
      }
1045 1047

	
1046 1048
      // Initialize arc maps
1047 1049
      for (int i = 0; i != _arc_num; ++i) {
1048 1050
        _flow[i] = 0;
1049 1051
        _state[i] = STATE_LOWER;
1050 1052
      }
1051 1053
      
1052 1054
      // Set data for the artificial root node
1053 1055
      _root = _node_num;
1054 1056
      _parent[_root] = -1;
1055 1057
      _pred[_root] = -1;
1056 1058
      _thread[_root] = 0;
1057 1059
      _rev_thread[0] = _root;
1058 1060
      _succ_num[_root] = _node_num + 1;
1059 1061
      _last_succ[_root] = _root - 1;
1060 1062
      _supply[_root] = -_sum_supply;
1061 1063
      _pi[_root] = 0;
1062 1064

	
1063 1065
      // Add artificial arcs and initialize the spanning tree data structure
1064 1066
      if (_sum_supply == 0) {
1065 1067
        // EQ supply constraints
1066 1068
        _search_arc_num = _arc_num;
1067 1069
        _all_arc_num = _arc_num + _node_num;
1068 1070
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1069 1071
          _parent[u] = _root;
1070 1072
          _pred[u] = e;
1071 1073
          _thread[u] = u + 1;
1072 1074
          _rev_thread[u + 1] = u;
1073 1075
          _succ_num[u] = 1;
1074 1076
          _last_succ[u] = u;
1075 1077
          _cap[e] = INF;
1076 1078
          _state[e] = STATE_TREE;
1077 1079
          if (_supply[u] >= 0) {
1078 1080
            _forward[u] = true;
1079 1081
            _pi[u] = 0;
1080 1082
            _source[e] = u;
1081 1083
            _target[e] = _root;
1082 1084
            _flow[e] = _supply[u];
1083 1085
            _cost[e] = 0;
1084 1086
          } else {
1085 1087
            _forward[u] = false;
1086 1088
            _pi[u] = ART_COST;
1087 1089
            _source[e] = _root;
1088 1090
            _target[e] = u;
1089 1091
            _flow[e] = -_supply[u];
1090 1092
            _cost[e] = ART_COST;
1091 1093
          }
1092 1094
        }
1093 1095
      }
1094 1096
      else if (_sum_supply > 0) {
1095 1097
        // LEQ supply constraints
1096 1098
        _search_arc_num = _arc_num + _node_num;
1097 1099
        int f = _arc_num + _node_num;
1098 1100
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1099 1101
          _parent[u] = _root;
1100 1102
          _thread[u] = u + 1;
1101 1103
          _rev_thread[u + 1] = u;
1102 1104
          _succ_num[u] = 1;
1103 1105
          _last_succ[u] = u;
1104 1106
          if (_supply[u] >= 0) {
1105 1107
            _forward[u] = true;
1106 1108
            _pi[u] = 0;
1107 1109
            _pred[u] = e;
1108 1110
            _source[e] = u;
1109 1111
            _target[e] = _root;
1110 1112
            _cap[e] = INF;
1111 1113
            _flow[e] = _supply[u];
1112 1114
            _cost[e] = 0;
1113 1115
            _state[e] = STATE_TREE;
1114 1116
          } else {
1115 1117
            _forward[u] = false;
1116 1118
            _pi[u] = ART_COST;
1117 1119
            _pred[u] = f;
1118 1120
            _source[f] = _root;
1119 1121
            _target[f] = u;
1120 1122
            _cap[f] = INF;
1121 1123
            _flow[f] = -_supply[u];
1122 1124
            _cost[f] = ART_COST;
1123 1125
            _state[f] = STATE_TREE;
1124 1126
            _source[e] = u;
1125 1127
            _target[e] = _root;
1126 1128
            _cap[e] = INF;
1127 1129
            _flow[e] = 0;
1128 1130
            _cost[e] = 0;
1129 1131
            _state[e] = STATE_LOWER;
1130 1132
            ++f;
1131 1133
          }
1132 1134
        }
1133 1135
        _all_arc_num = f;
1134 1136
      }
1135 1137
      else {
1136 1138
        // GEQ supply constraints
1137 1139
        _search_arc_num = _arc_num + _node_num;
1138 1140
        int f = _arc_num + _node_num;
1139 1141
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1140 1142
          _parent[u] = _root;
1141 1143
          _thread[u] = u + 1;
1142 1144
          _rev_thread[u + 1] = u;
1143 1145
          _succ_num[u] = 1;
1144 1146
          _last_succ[u] = u;
1145 1147
          if (_supply[u] <= 0) {
1146 1148
            _forward[u] = false;
1147 1149
            _pi[u] = 0;
1148 1150
            _pred[u] = e;
1149 1151
            _source[e] = _root;
1150 1152
            _target[e] = u;
1151 1153
            _cap[e] = INF;
1152 1154
            _flow[e] = -_supply[u];
1153 1155
            _cost[e] = 0;
1154 1156
            _state[e] = STATE_TREE;
1155 1157
          } else {
1156 1158
            _forward[u] = true;
1157 1159
            _pi[u] = -ART_COST;
1158 1160
            _pred[u] = f;
1159 1161
            _source[f] = u;
1160 1162
            _target[f] = _root;
1161 1163
            _cap[f] = INF;
1162 1164
            _flow[f] = _supply[u];
1163 1165
            _state[f] = STATE_TREE;
1164 1166
            _cost[f] = ART_COST;
1165 1167
            _source[e] = _root;
1166 1168
            _target[e] = u;
1167 1169
            _cap[e] = INF;
1168 1170
            _flow[e] = 0;
1169 1171
            _cost[e] = 0;
1170 1172
            _state[e] = STATE_LOWER;
1171 1173
            ++f;
1172 1174
          }
1173 1175
        }
1174 1176
        _all_arc_num = f;
1175 1177
      }
1176 1178

	
1177 1179
      return true;
1178 1180
    }
1179 1181

	
1180 1182
    // Find the join node
1181 1183
    void findJoinNode() {
1182 1184
      int u = _source[in_arc];
1183 1185
      int v = _target[in_arc];
1184 1186
      while (u != v) {
1185 1187
        if (_succ_num[u] < _succ_num[v]) {
1186 1188
          u = _parent[u];
1187 1189
        } else {
1188 1190
          v = _parent[v];
1189 1191
        }
1190 1192
      }
1191 1193
      join = u;
1192 1194
    }
1193 1195

	
1194 1196
    // Find the leaving arc of the cycle and returns true if the
1195 1197
    // leaving arc is not the same as the entering arc
1196 1198
    bool findLeavingArc() {
1197 1199
      // Initialize first and second nodes according to the direction
1198 1200
      // of the cycle
1199 1201
      if (_state[in_arc] == STATE_LOWER) {
1200 1202
        first  = _source[in_arc];
1201 1203
        second = _target[in_arc];
1202 1204
      } else {
1203 1205
        first  = _target[in_arc];
1204 1206
        second = _source[in_arc];
1205 1207
      }
1206 1208
      delta = _cap[in_arc];
1207 1209
      int result = 0;
1208 1210
      Value d;
1209 1211
      int e;
1210 1212

	
1211 1213
      // Search the cycle along the path form the first node to the root
1212 1214
      for (int u = first; u != join; u = _parent[u]) {
1213 1215
        e = _pred[u];
1214 1216
        d = _forward[u] ?
1215 1217
          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
1216 1218
        if (d < delta) {
1217 1219
          delta = d;
1218 1220
          u_out = u;
1219 1221
          result = 1;
1220 1222
        }
1221 1223
      }
1222 1224
      // Search the cycle along the path form the second node to the root
1223 1225
      for (int u = second; u != join; u = _parent[u]) {
1224 1226
        e = _pred[u];
1225 1227
        d = _forward[u] ? 
1226 1228
          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
1227 1229
        if (d <= delta) {
1228 1230
          delta = d;
1229 1231
          u_out = u;
1230 1232
          result = 2;
1231 1233
        }
1232 1234
      }
1233 1235

	
1234 1236
      if (result == 1) {
1235 1237
        u_in = first;
1236 1238
        v_in = second;
1237 1239
      } else {
1238 1240
        u_in = second;
1239 1241
        v_in = first;
1240 1242
      }
1241 1243
      return result != 0;
1242 1244
    }
1243 1245

	
1244 1246
    // Change _flow and _state vectors
1245 1247
    void changeFlow(bool change) {
1246 1248
      // Augment along the cycle
1247 1249
      if (delta > 0) {
1248 1250
        Value val = _state[in_arc] * delta;
1249 1251
        _flow[in_arc] += val;
1250 1252
        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
1251 1253
          _flow[_pred[u]] += _forward[u] ? -val : val;
1252 1254
        }
1253 1255
        for (int u = _target[in_arc]; u != join; u = _parent[u]) {
1254 1256
          _flow[_pred[u]] += _forward[u] ? val : -val;
1255 1257
        }
1256 1258
      }
1257 1259
      // Update the state of the entering and leaving arcs
1258 1260
      if (change) {
1259 1261
        _state[in_arc] = STATE_TREE;
1260 1262
        _state[_pred[u_out]] =
1261 1263
          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
1262 1264
      } else {
1263 1265
        _state[in_arc] = -_state[in_arc];
1264 1266
      }
1265 1267
    }
1266 1268

	
1267 1269
    // Update the tree structure
1268 1270
    void updateTreeStructure() {
1269 1271
      int u, w;
1270 1272
      int old_rev_thread = _rev_thread[u_out];
1271 1273
      int old_succ_num = _succ_num[u_out];
1272 1274
      int old_last_succ = _last_succ[u_out];
1273 1275
      v_out = _parent[u_out];
1274 1276

	
1275 1277
      u = _last_succ[u_in];  // the last successor of u_in
1276 1278
      right = _thread[u];    // the node after it
1277 1279

	
1278 1280
      // Handle the case when old_rev_thread equals to v_in
1279 1281
      // (it also means that join and v_out coincide)
1280 1282
      if (old_rev_thread == v_in) {
1281 1283
        last = _thread[_last_succ[u_out]];
1282 1284
      } else {
1283 1285
        last = _thread[v_in];
1284 1286
      }
1285 1287

	
1286 1288
      // Update _thread and _parent along the stem nodes (i.e. the nodes
1287 1289
      // between u_in and u_out, whose parent have to be changed)
1288 1290
      _thread[v_in] = stem = u_in;
1289 1291
      _dirty_revs.clear();
1290 1292
      _dirty_revs.push_back(v_in);
1291 1293
      par_stem = v_in;
1292 1294
      while (stem != u_out) {
1293 1295
        // Insert the next stem node into the thread list
1294 1296
        new_stem = _parent[stem];
1295 1297
        _thread[u] = new_stem;
1296 1298
        _dirty_revs.push_back(u);
1297 1299

	
1298 1300
        // Remove the subtree of stem from the thread list
1299 1301
        w = _rev_thread[stem];
1300 1302
        _thread[w] = right;
1301 1303
        _rev_thread[right] = w;
1302 1304

	
1303 1305
        // Change the parent node and shift stem nodes
1304 1306
        _parent[stem] = par_stem;
1305 1307
        par_stem = stem;
1306 1308
        stem = new_stem;
1307 1309

	
1308 1310
        // Update u and right
1309 1311
        u = _last_succ[stem] == _last_succ[par_stem] ?
1310 1312
          _rev_thread[par_stem] : _last_succ[stem];
1311 1313
        right = _thread[u];
1312 1314
      }
1313 1315
      _parent[u_out] = par_stem;
1314 1316
      _thread[u] = last;
1315 1317
      _rev_thread[last] = u;
1316 1318
      _last_succ[u_out] = u;
1317 1319

	
1318 1320
      // Remove the subtree of u_out from the thread list except for
1319 1321
      // the case when old_rev_thread equals to v_in
1320 1322
      // (it also means that join and v_out coincide)
1321 1323
      if (old_rev_thread != v_in) {
1322 1324
        _thread[old_rev_thread] = right;
1323 1325
        _rev_thread[right] = old_rev_thread;
1324 1326
      }
1325 1327

	
1326 1328
      // Update _rev_thread using the new _thread values
1327 1329
      for (int i = 0; i < int(_dirty_revs.size()); ++i) {
1328 1330
        u = _dirty_revs[i];
1329 1331
        _rev_thread[_thread[u]] = u;
1330 1332
      }
1331 1333

	
1332 1334
      // Update _pred, _forward, _last_succ and _succ_num for the
1333 1335
      // stem nodes from u_out to u_in
1334 1336
      int tmp_sc = 0, tmp_ls = _last_succ[u_out];
1335 1337
      u = u_out;
1336 1338
      while (u != u_in) {
1337 1339
        w = _parent[u];
1338 1340
        _pred[u] = _pred[w];
1339 1341
        _forward[u] = !_forward[w];
1340 1342
        tmp_sc += _succ_num[u] - _succ_num[w];
1341 1343
        _succ_num[u] = tmp_sc;
1342 1344
        _last_succ[w] = tmp_ls;
1343 1345
        u = w;
1344 1346
      }
1345 1347
      _pred[u_in] = in_arc;
1346 1348
      _forward[u_in] = (u_in == _source[in_arc]);
1347 1349
      _succ_num[u_in] = old_succ_num;
1348 1350

	
1349 1351
      // Set limits for updating _last_succ form v_in and v_out
1350 1352
      // towards the root
1351 1353
      int up_limit_in = -1;
1352 1354
      int up_limit_out = -1;
1353 1355
      if (_last_succ[join] == v_in) {
1354 1356
        up_limit_out = join;
1355 1357
      } else {
1356 1358
        up_limit_in = join;
1357 1359
      }
1358 1360

	
1359 1361
      // Update _last_succ from v_in towards the root
1360 1362
      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
1361 1363
           u = _parent[u]) {
1362 1364
        _last_succ[u] = _last_succ[u_out];
1363 1365
      }
1364 1366
      // Update _last_succ from v_out towards the root
1365 1367
      if (join != old_rev_thread && v_in != old_rev_thread) {
1366 1368
        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1367 1369
             u = _parent[u]) {
1368 1370
          _last_succ[u] = old_rev_thread;
1369 1371
        }
1370 1372
      } else {
1371 1373
        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1372 1374
             u = _parent[u]) {
1373 1375
          _last_succ[u] = _last_succ[u_out];
1374 1376
        }
1375 1377
      }
1376 1378

	
1377 1379
      // Update _succ_num from v_in to join
1378 1380
      for (u = v_in; u != join; u = _parent[u]) {
1379 1381
        _succ_num[u] += old_succ_num;
1380 1382
      }
1381 1383
      // Update _succ_num from v_out to join
1382 1384
      for (u = v_out; u != join; u = _parent[u]) {
1383 1385
        _succ_num[u] -= old_succ_num;
1384 1386
      }
1385 1387
    }
1386 1388

	
1387 1389
    // Update potentials
1388 1390
    void updatePotential() {
1389 1391
      Cost sigma = _forward[u_in] ?
1390 1392
        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1391 1393
        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1392 1394
      // Update potentials in the subtree, which has been moved
1393 1395
      int end = _thread[_last_succ[u_in]];
1394 1396
      for (int u = u_in; u != end; u = _thread[u]) {
1395 1397
        _pi[u] += sigma;
1396 1398
      }
1397 1399
    }
1398 1400

	
1399 1401
    // Execute the algorithm
1400 1402
    ProblemType start(PivotRule pivot_rule) {
1401 1403
      // Select the pivot rule implementation
1402 1404
      switch (pivot_rule) {
1403 1405
        case FIRST_ELIGIBLE:
1404 1406
          return start<FirstEligiblePivotRule>();
1405 1407
        case BEST_ELIGIBLE:
1406 1408
          return start<BestEligiblePivotRule>();
1407 1409
        case BLOCK_SEARCH:
1408 1410
          return start<BlockSearchPivotRule>();
1409 1411
        case CANDIDATE_LIST:
1410 1412
          return start<CandidateListPivotRule>();
1411 1413
        case ALTERING_LIST:
1412 1414
          return start<AlteringListPivotRule>();
1413 1415
      }
1414 1416
      return INFEASIBLE; // avoid warning
1415 1417
    }
1416 1418

	
1417 1419
    template <typename PivotRuleImpl>
1418 1420
    ProblemType start() {
1419 1421
      PivotRuleImpl pivot(*this);
1420 1422

	
1421 1423
      // Execute the Network Simplex algorithm
1422 1424
      while (pivot.findEnteringArc()) {
1423 1425
        findJoinNode();
1424 1426
        bool change = findLeavingArc();
1425 1427
        if (delta >= INF) return UNBOUNDED;
1426 1428
        changeFlow(change);
1427 1429
        if (change) {
1428 1430
          updateTreeStructure();
1429 1431
          updatePotential();
1430 1432
        }
1431 1433
      }
1432 1434
      
1433 1435
      // Check feasibility
1434 1436
      for (int e = _search_arc_num; e != _all_arc_num; ++e) {
1435 1437
        if (_flow[e] != 0) return INFEASIBLE;
1436 1438
      }
1437 1439

	
1438 1440
      // Transform the solution and the supply map to the original form
1439 1441
      if (_have_lower) {
1440 1442
        for (int i = 0; i != _arc_num; ++i) {
1441 1443
          Value c = _lower[i];
1442 1444
          if (c != 0) {
1443 1445
            _flow[i] += c;
1444 1446
            _supply[_source[i]] += c;
1445 1447
            _supply[_target[i]] -= c;
1446 1448
          }
1447 1449
        }
1448 1450
      }
1449 1451
      
1450 1452
      // Shift potentials to meet the requirements of the GEQ/LEQ type
1451 1453
      // optimality conditions
1452 1454
      if (_sum_supply == 0) {
1453 1455
        if (_stype == GEQ) {
1454 1456
          Cost max_pot = std::numeric_limits<Cost>::min();
1455 1457
          for (int i = 0; i != _node_num; ++i) {
1456 1458
            if (_pi[i] > max_pot) max_pot = _pi[i];
1457 1459
          }
1458 1460
          if (max_pot > 0) {
1459 1461
            for (int i = 0; i != _node_num; ++i)
1460 1462
              _pi[i] -= max_pot;
1461 1463
          }
1462 1464
        } else {
1463 1465
          Cost min_pot = std::numeric_limits<Cost>::max();
1464 1466
          for (int i = 0; i != _node_num; ++i) {
1465 1467
            if (_pi[i] < min_pot) min_pot = _pi[i];
1466 1468
          }
1467 1469
          if (min_pot < 0) {
1468 1470
            for (int i = 0; i != _node_num; ++i)
1469 1471
              _pi[i] -= min_pot;
1470 1472
          }
1471 1473
        }
1472 1474
      }
1473 1475

	
1474 1476
      return OPTIMAL;
1475 1477
    }
1476 1478

	
1477 1479
  }; //class NetworkSimplex
1478 1480

	
1479 1481
  ///@}
1480 1482

	
1481 1483
} //namespace lemon
1482 1484

	
1483 1485
#endif //LEMON_NETWORK_SIMPLEX_H
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_PREFLOW_H
20 20
#define LEMON_PREFLOW_H
21 21

	
22 22
#include <lemon/tolerance.h>
23 23
#include <lemon/elevator.h>
24 24

	
25 25
/// \file
26 26
/// \ingroup max_flow
27 27
/// \brief Implementation of the preflow algorithm.
28 28

	
29 29
namespace lemon {
30 30

	
31 31
  /// \brief Default traits class of Preflow class.
32 32
  ///
33 33
  /// Default traits class of Preflow class.
34 34
  /// \tparam GR Digraph type.
35 35
  /// \tparam CAP Capacity map type.
36 36
  template <typename GR, typename CAP>
37 37
  struct PreflowDefaultTraits {
38 38

	
39 39
    /// \brief The type of the digraph the algorithm runs on.
40 40
    typedef GR Digraph;
41 41

	
42 42
    /// \brief The type of the map that stores the arc capacities.
43 43
    ///
44 44
    /// The type of the map that stores the arc capacities.
45 45
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
46 46
    typedef CAP CapacityMap;
47 47

	
48 48
    /// \brief The type of the flow values.
49 49
    typedef typename CapacityMap::Value Value;
50 50

	
51 51
    /// \brief The type of the map that stores the flow values.
52 52
    ///
53 53
    /// The type of the map that stores the flow values.
54 54
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
55 55
#ifdef DOXYGEN
56 56
    typedef GR::ArcMap<Value> FlowMap;
57 57
#else
58 58
    typedef typename Digraph::template ArcMap<Value> FlowMap;
59 59
#endif
60 60

	
61 61
    /// \brief Instantiates a FlowMap.
62 62
    ///
63 63
    /// This function instantiates a \ref FlowMap.
64 64
    /// \param digraph The digraph for which we would like to define
65 65
    /// the flow map.
66 66
    static FlowMap* createFlowMap(const Digraph& digraph) {
67 67
      return new FlowMap(digraph);
68 68
    }
69 69

	
70 70
    /// \brief The elevator type used by Preflow algorithm.
71 71
    ///
72 72
    /// The elevator type used by Preflow algorithm.
73 73
    ///
74 74
    /// \sa Elevator, LinkedElevator
75 75
#ifdef DOXYGEN
76 76
    typedef lemon::Elevator<GR, GR::Node> Elevator;
77 77
#else
78 78
    typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
79 79
#endif
80 80

	
81 81
    /// \brief Instantiates an Elevator.
82 82
    ///
83 83
    /// This function instantiates an \ref Elevator.
84 84
    /// \param digraph The digraph for which we would like to define
85 85
    /// the elevator.
86 86
    /// \param max_level The maximum level of the elevator.
87 87
    static Elevator* createElevator(const Digraph& digraph, int max_level) {
88 88
      return new Elevator(digraph, max_level);
89 89
    }
90 90

	
91 91
    /// \brief The tolerance used by the algorithm
92 92
    ///
93 93
    /// The tolerance used by the algorithm to handle inexact computation.
94 94
    typedef lemon::Tolerance<Value> Tolerance;
95 95

	
96 96
  };
97 97

	
98 98

	
99 99
  /// \ingroup max_flow
100 100
  ///
101 101
  /// \brief %Preflow algorithm class.
102 102
  ///
103 103
  /// This class provides an implementation of Goldberg-Tarjan's \e preflow
104 104
  /// \e push-relabel algorithm producing a \ref max_flow
105
  /// "flow of maximum value" in a digraph.
105
  /// "flow of maximum value" in a digraph \ref clrs01algorithms,
106
  /// \ref amo93networkflows, \ref goldberg88newapproach.
106 107
  /// The preflow algorithms are the fastest known maximum
107 108
  /// flow algorithms. The current implementation uses a mixture of the
108 109
  /// \e "highest label" and the \e "bound decrease" heuristics.
109 110
  /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
110 111
  ///
111 112
  /// The algorithm consists of two phases. After the first phase
112 113
  /// the maximum flow value and the minimum cut is obtained. The
113 114
  /// second phase constructs a feasible maximum flow on each arc.
114 115
  ///
115 116
  /// \tparam GR The type of the digraph the algorithm runs on.
116 117
  /// \tparam CAP The type of the capacity map. The default map
117 118
  /// type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
118 119
#ifdef DOXYGEN
119 120
  template <typename GR, typename CAP, typename TR>
120 121
#else
121 122
  template <typename GR,
122 123
            typename CAP = typename GR::template ArcMap<int>,
123 124
            typename TR = PreflowDefaultTraits<GR, CAP> >
124 125
#endif
125 126
  class Preflow {
126 127
  public:
127 128

	
128 129
    ///The \ref PreflowDefaultTraits "traits class" of the algorithm.
129 130
    typedef TR Traits;
130 131
    ///The type of the digraph the algorithm runs on.
131 132
    typedef typename Traits::Digraph Digraph;
132 133
    ///The type of the capacity map.
133 134
    typedef typename Traits::CapacityMap CapacityMap;
134 135
    ///The type of the flow values.
135 136
    typedef typename Traits::Value Value;
136 137

	
137 138
    ///The type of the flow map.
138 139
    typedef typename Traits::FlowMap FlowMap;
139 140
    ///The type of the elevator.
140 141
    typedef typename Traits::Elevator Elevator;
141 142
    ///The type of the tolerance.
142 143
    typedef typename Traits::Tolerance Tolerance;
143 144

	
144 145
  private:
145 146

	
146 147
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
147 148

	
148 149
    const Digraph& _graph;
149 150
    const CapacityMap* _capacity;
150 151

	
151 152
    int _node_num;
152 153

	
153 154
    Node _source, _target;
154 155

	
155 156
    FlowMap* _flow;
156 157
    bool _local_flow;
157 158

	
158 159
    Elevator* _level;
159 160
    bool _local_level;
160 161

	
161 162
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
162 163
    ExcessMap* _excess;
163 164

	
164 165
    Tolerance _tolerance;
165 166

	
166 167
    bool _phase;
167 168

	
168 169

	
169 170
    void createStructures() {
170 171
      _node_num = countNodes(_graph);
171 172

	
172 173
      if (!_flow) {
173 174
        _flow = Traits::createFlowMap(_graph);
174 175
        _local_flow = true;
175 176
      }
176 177
      if (!_level) {
177 178
        _level = Traits::createElevator(_graph, _node_num);
178 179
        _local_level = true;
179 180
      }
180 181
      if (!_excess) {
181 182
        _excess = new ExcessMap(_graph);
182 183
      }
183 184
    }
184 185

	
185 186
    void destroyStructures() {
186 187
      if (_local_flow) {
187 188
        delete _flow;
188 189
      }
189 190
      if (_local_level) {
190 191
        delete _level;
191 192
      }
192 193
      if (_excess) {
193 194
        delete _excess;
194 195
      }
195 196
    }
196 197

	
197 198
  public:
198 199

	
199 200
    typedef Preflow Create;
200 201

	
201 202
    ///\name Named Template Parameters
202 203

	
203 204
    ///@{
204 205

	
205 206
    template <typename T>
206 207
    struct SetFlowMapTraits : public Traits {
207 208
      typedef T FlowMap;
208 209
      static FlowMap *createFlowMap(const Digraph&) {
209 210
        LEMON_ASSERT(false, "FlowMap is not initialized");
210 211
        return 0; // ignore warnings
211 212
      }
212 213
    };
213 214

	
214 215
    /// \brief \ref named-templ-param "Named parameter" for setting
215 216
    /// FlowMap type
216 217
    ///
217 218
    /// \ref named-templ-param "Named parameter" for setting FlowMap
218 219
    /// type.
219 220
    template <typename T>
220 221
    struct SetFlowMap
221 222
      : public Preflow<Digraph, CapacityMap, SetFlowMapTraits<T> > {
222 223
      typedef Preflow<Digraph, CapacityMap,
223 224
                      SetFlowMapTraits<T> > Create;
224 225
    };
225 226

	
226 227
    template <typename T>
227 228
    struct SetElevatorTraits : public Traits {
228 229
      typedef T Elevator;
229 230
      static Elevator *createElevator(const Digraph&, int) {
230 231
        LEMON_ASSERT(false, "Elevator is not initialized");
231 232
        return 0; // ignore warnings
232 233
      }
233 234
    };
234 235

	
235 236
    /// \brief \ref named-templ-param "Named parameter" for setting
236 237
    /// Elevator type
237 238
    ///
238 239
    /// \ref named-templ-param "Named parameter" for setting Elevator
239 240
    /// type. If this named parameter is used, then an external
240 241
    /// elevator object must be passed to the algorithm using the
241 242
    /// \ref elevator(Elevator&) "elevator()" function before calling
242 243
    /// \ref run() or \ref init().
243 244
    /// \sa SetStandardElevator
244 245
    template <typename T>
245 246
    struct SetElevator
246 247
      : public Preflow<Digraph, CapacityMap, SetElevatorTraits<T> > {
247 248
      typedef Preflow<Digraph, CapacityMap,
248 249
                      SetElevatorTraits<T> > Create;
249 250
    };
250 251

	
251 252
    template <typename T>
252 253
    struct SetStandardElevatorTraits : public Traits {
253 254
      typedef T Elevator;
254 255
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
255 256
        return new Elevator(digraph, max_level);
256 257
      }
257 258
    };
258 259

	
259 260
    /// \brief \ref named-templ-param "Named parameter" for setting
260 261
    /// Elevator type with automatic allocation
261 262
    ///
262 263
    /// \ref named-templ-param "Named parameter" for setting Elevator
263 264
    /// type with automatic allocation.
264 265
    /// The Elevator should have standard constructor interface to be
265 266
    /// able to automatically created by the algorithm (i.e. the
266 267
    /// digraph and the maximum level should be passed to it).
267 268
    /// However an external elevator object could also be passed to the
268 269
    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
269 270
    /// before calling \ref run() or \ref init().
270 271
    /// \sa SetElevator
271 272
    template <typename T>
272 273
    struct SetStandardElevator
273 274
      : public Preflow<Digraph, CapacityMap,
274 275
                       SetStandardElevatorTraits<T> > {
275 276
      typedef Preflow<Digraph, CapacityMap,
276 277
                      SetStandardElevatorTraits<T> > Create;
277 278
    };
278 279

	
279 280
    /// @}
280 281

	
281 282
  protected:
282 283

	
283 284
    Preflow() {}
284 285

	
285 286
  public:
286 287

	
287 288

	
288 289
    /// \brief The constructor of the class.
289 290
    ///
290 291
    /// The constructor of the class.
291 292
    /// \param digraph The digraph the algorithm runs on.
292 293
    /// \param capacity The capacity of the arcs.
293 294
    /// \param source The source node.
294 295
    /// \param target The target node.
295 296
    Preflow(const Digraph& digraph, const CapacityMap& capacity,
296 297
            Node source, Node target)
297 298
      : _graph(digraph), _capacity(&capacity),
298 299
        _node_num(0), _source(source), _target(target),
299 300
        _flow(0), _local_flow(false),
300 301
        _level(0), _local_level(false),
301 302
        _excess(0), _tolerance(), _phase() {}
302 303

	
303 304
    /// \brief Destructor.
304 305
    ///
305 306
    /// Destructor.
306 307
    ~Preflow() {
307 308
      destroyStructures();
308 309
    }
309 310

	
310 311
    /// \brief Sets the capacity map.
311 312
    ///
312 313
    /// Sets the capacity map.
313 314
    /// \return <tt>(*this)</tt>
314 315
    Preflow& capacityMap(const CapacityMap& map) {
315 316
      _capacity = &map;
316 317
      return *this;
317 318
    }
318 319

	
319 320
    /// \brief Sets the flow map.
320 321
    ///
321 322
    /// Sets the flow map.
322 323
    /// If you don't use this function before calling \ref run() or
323 324
    /// \ref init(), an instance will be allocated automatically.
324 325
    /// The destructor deallocates this automatically allocated map,
325 326
    /// of course.
326 327
    /// \return <tt>(*this)</tt>
327 328
    Preflow& flowMap(FlowMap& map) {
328 329
      if (_local_flow) {
329 330
        delete _flow;
330 331
        _local_flow = false;
331 332
      }
332 333
      _flow = &map;
333 334
      return *this;
334 335
    }
335 336

	
336 337
    /// \brief Sets the source node.
337 338
    ///
338 339
    /// Sets the source node.
339 340
    /// \return <tt>(*this)</tt>
340 341
    Preflow& source(const Node& node) {
341 342
      _source = node;
342 343
      return *this;
343 344
    }
344 345

	
345 346
    /// \brief Sets the target node.
346 347
    ///
347 348
    /// Sets the target node.
348 349
    /// \return <tt>(*this)</tt>
349 350
    Preflow& target(const Node& node) {
350 351
      _target = node;
351 352
      return *this;
352 353
    }
353 354

	
354 355
    /// \brief Sets the elevator used by algorithm.
355 356
    ///
356 357
    /// Sets the elevator used by algorithm.
357 358
    /// If you don't use this function before calling \ref run() or
358 359
    /// \ref init(), an instance will be allocated automatically.
359 360
    /// The destructor deallocates this automatically allocated elevator,
360 361
    /// of course.
361 362
    /// \return <tt>(*this)</tt>
362 363
    Preflow& elevator(Elevator& elevator) {
363 364
      if (_local_level) {
364 365
        delete _level;
365 366
        _local_level = false;
366 367
      }
367 368
      _level = &elevator;
368 369
      return *this;
369 370
    }
370 371

	
371 372
    /// \brief Returns a const reference to the elevator.
372 373
    ///
373 374
    /// Returns a const reference to the elevator.
374 375
    ///
375 376
    /// \pre Either \ref run() or \ref init() must be called before
376 377
    /// using this function.
377 378
    const Elevator& elevator() const {
378 379
      return *_level;
379 380
    }
380 381

	
381 382
    /// \brief Sets the tolerance used by the algorithm.
382 383
    ///
383 384
    /// Sets the tolerance object used by the algorithm.
384 385
    /// \return <tt>(*this)</tt>
385 386
    Preflow& tolerance(const Tolerance& tolerance) {
386 387
      _tolerance = tolerance;
387 388
      return *this;
388 389
    }
389 390

	
390 391
    /// \brief Returns a const reference to the tolerance.
391 392
    ///
392 393
    /// Returns a const reference to the tolerance object used by
393 394
    /// the algorithm.
394 395
    const Tolerance& tolerance() const {
395 396
      return _tolerance;
396 397
    }
397 398

	
398 399
    /// \name Execution Control
399 400
    /// The simplest way to execute the preflow algorithm is to use
400 401
    /// \ref run() or \ref runMinCut().\n
401 402
    /// If you need better control on the initial solution or the execution,
402 403
    /// you have to call one of the \ref init() functions first, then
403 404
    /// \ref startFirstPhase() and if you need it \ref startSecondPhase().
404 405

	
405 406
    ///@{
406 407

	
407 408
    /// \brief Initializes the internal data structures.
408 409
    ///
409 410
    /// Initializes the internal data structures and sets the initial
410 411
    /// flow to zero on each arc.
411 412
    void init() {
412 413
      createStructures();
413 414

	
414 415
      _phase = true;
415 416
      for (NodeIt n(_graph); n != INVALID; ++n) {
416 417
        (*_excess)[n] = 0;
417 418
      }
418 419

	
419 420
      for (ArcIt e(_graph); e != INVALID; ++e) {
420 421
        _flow->set(e, 0);
421 422
      }
422 423

	
423 424
      typename Digraph::template NodeMap<bool> reached(_graph, false);
424 425

	
425 426
      _level->initStart();
426 427
      _level->initAddItem(_target);
427 428

	
428 429
      std::vector<Node> queue;
429 430
      reached[_source] = true;
430 431

	
431 432
      queue.push_back(_target);
432 433
      reached[_target] = true;
433 434
      while (!queue.empty()) {
434 435
        _level->initNewLevel();
435 436
        std::vector<Node> nqueue;
436 437
        for (int i = 0; i < int(queue.size()); ++i) {
437 438
          Node n = queue[i];
438 439
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
439 440
            Node u = _graph.source(e);
440 441
            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
441 442
              reached[u] = true;
442 443
              _level->initAddItem(u);
443 444
              nqueue.push_back(u);
444 445
            }
445 446
          }
446 447
        }
447 448
        queue.swap(nqueue);
448 449
      }
449 450
      _level->initFinish();
450 451

	
451 452
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
452 453
        if (_tolerance.positive((*_capacity)[e])) {
453 454
          Node u = _graph.target(e);
454 455
          if ((*_level)[u] == _level->maxLevel()) continue;
455 456
          _flow->set(e, (*_capacity)[e]);
456 457
          (*_excess)[u] += (*_capacity)[e];
457 458
          if (u != _target && !_level->active(u)) {
458 459
            _level->activate(u);
459 460
          }
460 461
        }
461 462
      }
462 463
    }
463 464

	
464 465
    /// \brief Initializes the internal data structures using the
465 466
    /// given flow map.
466 467
    ///
467 468
    /// Initializes the internal data structures and sets the initial
468 469
    /// flow to the given \c flowMap. The \c flowMap should contain a
469 470
    /// flow or at least a preflow, i.e. at each node excluding the
470 471
    /// source node the incoming flow should greater or equal to the
471 472
    /// outgoing flow.
472 473
    /// \return \c false if the given \c flowMap is not a preflow.
473 474
    template <typename FlowMap>
474 475
    bool init(const FlowMap& flowMap) {
475 476
      createStructures();
476 477

	
477 478
      for (ArcIt e(_graph); e != INVALID; ++e) {
478 479
        _flow->set(e, flowMap[e]);
479 480
      }
480 481

	
481 482
      for (NodeIt n(_graph); n != INVALID; ++n) {
482 483
        Value excess = 0;
483 484
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
484 485
          excess += (*_flow)[e];
485 486
        }
486 487
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
487 488
          excess -= (*_flow)[e];
488 489
        }
489 490
        if (excess < 0 && n != _source) return false;
490 491
        (*_excess)[n] = excess;
491 492
      }
492 493

	
493 494
      typename Digraph::template NodeMap<bool> reached(_graph, false);
494 495

	
495 496
      _level->initStart();
496 497
      _level->initAddItem(_target);
497 498

	
498 499
      std::vector<Node> queue;
499 500
      reached[_source] = true;
500 501

	
501 502
      queue.push_back(_target);
502 503
      reached[_target] = true;
503 504
      while (!queue.empty()) {
504 505
        _level->initNewLevel();
505 506
        std::vector<Node> nqueue;
506 507
        for (int i = 0; i < int(queue.size()); ++i) {
507 508
          Node n = queue[i];
508 509
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
509 510
            Node u = _graph.source(e);
510 511
            if (!reached[u] &&
511 512
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
512 513
              reached[u] = true;
513 514
              _level->initAddItem(u);
514 515
              nqueue.push_back(u);
515 516
            }
516 517
          }
517 518
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
518 519
            Node v = _graph.target(e);
519 520
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
520 521
              reached[v] = true;
521 522
              _level->initAddItem(v);
522 523
              nqueue.push_back(v);
523 524
            }
524 525
          }
525 526
        }
526 527
        queue.swap(nqueue);
527 528
      }
528 529
      _level->initFinish();
529 530

	
530 531
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
531 532
        Value rem = (*_capacity)[e] - (*_flow)[e];
532 533
        if (_tolerance.positive(rem)) {
533 534
          Node u = _graph.target(e);
534 535
          if ((*_level)[u] == _level->maxLevel()) continue;
535 536
          _flow->set(e, (*_capacity)[e]);
536 537
          (*_excess)[u] += rem;
537 538
          if (u != _target && !_level->active(u)) {
538 539
            _level->activate(u);
539 540
          }
540 541
        }
541 542
      }
542 543
      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
543 544
        Value rem = (*_flow)[e];
544 545
        if (_tolerance.positive(rem)) {
545 546
          Node v = _graph.source(e);
546 547
          if ((*_level)[v] == _level->maxLevel()) continue;
547 548
          _flow->set(e, 0);
548 549
          (*_excess)[v] += rem;
549 550
          if (v != _target && !_level->active(v)) {
550 551
            _level->activate(v);
551 552
          }
552 553
        }
553 554
      }
554 555
      return true;
555 556
    }
556 557

	
557 558
    /// \brief Starts the first phase of the preflow algorithm.
558 559
    ///
559 560
    /// The preflow algorithm consists of two phases, this method runs
560 561
    /// the first phase. After the first phase the maximum flow value
561 562
    /// and a minimum value cut can already be computed, although a
562 563
    /// maximum flow is not yet obtained. So after calling this method
563 564
    /// \ref flowValue() returns the value of a maximum flow and \ref
564 565
    /// minCut() returns a minimum cut.
565 566
    /// \pre One of the \ref init() functions must be called before
566 567
    /// using this function.
567 568
    void startFirstPhase() {
568 569
      _phase = true;
569 570

	
570 571
      Node n = _level->highestActive();
571 572
      int level = _level->highestActiveLevel();
572 573
      while (n != INVALID) {
573 574
        int num = _node_num;
574 575

	
575 576
        while (num > 0 && n != INVALID) {
576 577
          Value excess = (*_excess)[n];
577 578
          int new_level = _level->maxLevel();
578 579

	
579 580
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
580 581
            Value rem = (*_capacity)[e] - (*_flow)[e];
581 582
            if (!_tolerance.positive(rem)) continue;
582 583
            Node v = _graph.target(e);
583 584
            if ((*_level)[v] < level) {
584 585
              if (!_level->active(v) && v != _target) {
585 586
                _level->activate(v);
586 587
              }
587 588
              if (!_tolerance.less(rem, excess)) {
588 589
                _flow->set(e, (*_flow)[e] + excess);
589 590
                (*_excess)[v] += excess;
590 591
                excess = 0;
591 592
                goto no_more_push_1;
592 593
              } else {
593 594
                excess -= rem;
594 595
                (*_excess)[v] += rem;
595 596
                _flow->set(e, (*_capacity)[e]);
596 597
              }
597 598
            } else if (new_level > (*_level)[v]) {
598 599
              new_level = (*_level)[v];
599 600
            }
600 601
          }
601 602

	
602 603
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
603 604
            Value rem = (*_flow)[e];
604 605
            if (!_tolerance.positive(rem)) continue;
605 606
            Node v = _graph.source(e);
606 607
            if ((*_level)[v] < level) {
607 608
              if (!_level->active(v) && v != _target) {
608 609
                _level->activate(v);
609 610
              }
610 611
              if (!_tolerance.less(rem, excess)) {
611 612
                _flow->set(e, (*_flow)[e] - excess);
612 613
                (*_excess)[v] += excess;
613 614
                excess = 0;
614 615
                goto no_more_push_1;
615 616
              } else {
616 617
                excess -= rem;
617 618
                (*_excess)[v] += rem;
618 619
                _flow->set(e, 0);
619 620
              }
620 621
            } else if (new_level > (*_level)[v]) {
621 622
              new_level = (*_level)[v];
622 623
            }
623 624
          }
624 625

	
625 626
        no_more_push_1:
626 627

	
627 628
          (*_excess)[n] = excess;
628 629

	
629 630
          if (excess != 0) {
630 631
            if (new_level + 1 < _level->maxLevel()) {
631 632
              _level->liftHighestActive(new_level + 1);
632 633
            } else {
633 634
              _level->liftHighestActiveToTop();
634 635
            }
635 636
            if (_level->emptyLevel(level)) {
636 637
              _level->liftToTop(level);
637 638
            }
638 639
          } else {
639 640
            _level->deactivate(n);
640 641
          }
641 642

	
642 643
          n = _level->highestActive();
643 644
          level = _level->highestActiveLevel();
644 645
          --num;
645 646
        }
646 647

	
647 648
        num = _node_num * 20;
648 649
        while (num > 0 && n != INVALID) {
649 650
          Value excess = (*_excess)[n];
650 651
          int new_level = _level->maxLevel();
651 652

	
652 653
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
653 654
            Value rem = (*_capacity)[e] - (*_flow)[e];
654 655
            if (!_tolerance.positive(rem)) continue;
655 656
            Node v = _graph.target(e);
656 657
            if ((*_level)[v] < level) {
657 658
              if (!_level->active(v) && v != _target) {
658 659
                _level->activate(v);
659 660
              }
660 661
              if (!_tolerance.less(rem, excess)) {
661 662
                _flow->set(e, (*_flow)[e] + excess);
662 663
                (*_excess)[v] += excess;
663 664
                excess = 0;
664 665
                goto no_more_push_2;
665 666
              } else {
666 667
                excess -= rem;
667 668
                (*_excess)[v] += rem;
668 669
                _flow->set(e, (*_capacity)[e]);
669 670
              }
670 671
            } else if (new_level > (*_level)[v]) {
671 672
              new_level = (*_level)[v];
672 673
            }
673 674
          }
674 675

	
675 676
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
676 677
            Value rem = (*_flow)[e];
677 678
            if (!_tolerance.positive(rem)) continue;
678 679
            Node v = _graph.source(e);
679 680
            if ((*_level)[v] < level) {
680 681
              if (!_level->active(v) && v != _target) {
681 682
                _level->activate(v);
682 683
              }
683 684
              if (!_tolerance.less(rem, excess)) {
684 685
                _flow->set(e, (*_flow)[e] - excess);
685 686
                (*_excess)[v] += excess;
686 687
                excess = 0;
687 688
                goto no_more_push_2;
688 689
              } else {
689 690
                excess -= rem;
690 691
                (*_excess)[v] += rem;
691 692
                _flow->set(e, 0);
692 693
              }
693 694
            } else if (new_level > (*_level)[v]) {
694 695
              new_level = (*_level)[v];
695 696
            }
696 697
          }
697 698

	
698 699
        no_more_push_2:
699 700

	
700 701
          (*_excess)[n] = excess;
701 702

	
702 703
          if (excess != 0) {
703 704
            if (new_level + 1 < _level->maxLevel()) {
704 705
              _level->liftActiveOn(level, new_level + 1);
705 706
            } else {
706 707
              _level->liftActiveToTop(level);
707 708
            }
708 709
            if (_level->emptyLevel(level)) {
709 710
              _level->liftToTop(level);
710 711
            }
711 712
          } else {
712 713
            _level->deactivate(n);
713 714
          }
714 715

	
715 716
          while (level >= 0 && _level->activeFree(level)) {
716 717
            --level;
717 718
          }
718 719
          if (level == -1) {
719 720
            n = _level->highestActive();
720 721
            level = _level->highestActiveLevel();
721 722
          } else {
722 723
            n = _level->activeOn(level);
723 724
          }
724 725
          --num;
725 726
        }
726 727
      }
727 728
    }
728 729

	
729 730
    /// \brief Starts the second phase of the preflow algorithm.
730 731
    ///
731 732
    /// The preflow algorithm consists of two phases, this method runs
732 733
    /// the second phase. After calling one of the \ref init() functions
733 734
    /// and \ref startFirstPhase() and then \ref startSecondPhase(),
734 735
    /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
735 736
    /// value of a maximum flow, \ref minCut() returns a minimum cut
736 737
    /// \pre One of the \ref init() functions and \ref startFirstPhase()
737 738
    /// must be called before using this function.
738 739
    void startSecondPhase() {
739 740
      _phase = false;
740 741

	
741 742
      typename Digraph::template NodeMap<bool> reached(_graph);
742 743
      for (NodeIt n(_graph); n != INVALID; ++n) {
743 744
        reached[n] = (*_level)[n] < _level->maxLevel();
744 745
      }
745 746

	
746 747
      _level->initStart();
747 748
      _level->initAddItem(_source);
748 749

	
749 750
      std::vector<Node> queue;
750 751
      queue.push_back(_source);
751 752
      reached[_source] = true;
752 753

	
753 754
      while (!queue.empty()) {
754 755
        _level->initNewLevel();
755 756
        std::vector<Node> nqueue;
756 757
        for (int i = 0; i < int(queue.size()); ++i) {
757 758
          Node n = queue[i];
758 759
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
759 760
            Node v = _graph.target(e);
760 761
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
761 762
              reached[v] = true;
762 763
              _level->initAddItem(v);
763 764
              nqueue.push_back(v);
764 765
            }
765 766
          }
766 767
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
767 768
            Node u = _graph.source(e);
768 769
            if (!reached[u] &&
769 770
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
770 771
              reached[u] = true;
771 772
              _level->initAddItem(u);
772 773
              nqueue.push_back(u);
773 774
            }
774 775
          }
775 776
        }
776 777
        queue.swap(nqueue);
777 778
      }
778 779
      _level->initFinish();
779 780

	
780 781
      for (NodeIt n(_graph); n != INVALID; ++n) {
781 782
        if (!reached[n]) {
782 783
          _level->dirtyTopButOne(n);
783 784
        } else if ((*_excess)[n] > 0 && _target != n) {
784 785
          _level->activate(n);
785 786
        }
786 787
      }
787 788

	
788 789
      Node n;
789 790
      while ((n = _level->highestActive()) != INVALID) {
790 791
        Value excess = (*_excess)[n];
791 792
        int level = _level->highestActiveLevel();
792 793
        int new_level = _level->maxLevel();
793 794

	
794 795
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
795 796
          Value rem = (*_capacity)[e] - (*_flow)[e];
796 797
          if (!_tolerance.positive(rem)) continue;
797 798
          Node v = _graph.target(e);
798 799
          if ((*_level)[v] < level) {
799 800
            if (!_level->active(v) && v != _source) {
800 801
              _level->activate(v);
801 802
            }
802 803
            if (!_tolerance.less(rem, excess)) {
803 804
              _flow->set(e, (*_flow)[e] + excess);
804 805
              (*_excess)[v] += excess;
805 806
              excess = 0;
806 807
              goto no_more_push;
807 808
            } else {
808 809
              excess -= rem;
809 810
              (*_excess)[v] += rem;
810 811
              _flow->set(e, (*_capacity)[e]);
811 812
            }
812 813
          } else if (new_level > (*_level)[v]) {
813 814
            new_level = (*_level)[v];
814 815
          }
815 816
        }
816 817

	
817 818
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
818 819
          Value rem = (*_flow)[e];
819 820
          if (!_tolerance.positive(rem)) continue;
820 821
          Node v = _graph.source(e);
821 822
          if ((*_level)[v] < level) {
822 823
            if (!_level->active(v) && v != _source) {
823 824
              _level->activate(v);
824 825
            }
825 826
            if (!_tolerance.less(rem, excess)) {
826 827
              _flow->set(e, (*_flow)[e] - excess);
827 828
              (*_excess)[v] += excess;
828 829
              excess = 0;
829 830
              goto no_more_push;
830 831
            } else {
831 832
              excess -= rem;
832 833
              (*_excess)[v] += rem;
833 834
              _flow->set(e, 0);
834 835
            }
835 836
          } else if (new_level > (*_level)[v]) {
836 837
            new_level = (*_level)[v];
837 838
          }
838 839
        }
839 840

	
840 841
      no_more_push:
841 842

	
842 843
        (*_excess)[n] = excess;
843 844

	
844 845
        if (excess != 0) {
845 846
          if (new_level + 1 < _level->maxLevel()) {
846 847
            _level->liftHighestActive(new_level + 1);
847 848
          } else {
848 849
            // Calculation error
849 850
            _level->liftHighestActiveToTop();
850 851
          }
851 852
          if (_level->emptyLevel(level)) {
852 853
            // Calculation error
853 854
            _level->liftToTop(level);
854 855
          }
855 856
        } else {
856 857
          _level->deactivate(n);
857 858
        }
858 859

	
859 860
      }
860 861
    }
861 862

	
862 863
    /// \brief Runs the preflow algorithm.
863 864
    ///
864 865
    /// Runs the preflow algorithm.
865 866
    /// \note pf.run() is just a shortcut of the following code.
866 867
    /// \code
867 868
    ///   pf.init();
868 869
    ///   pf.startFirstPhase();
869 870
    ///   pf.startSecondPhase();
870 871
    /// \endcode
871 872
    void run() {
872 873
      init();
873 874
      startFirstPhase();
874 875
      startSecondPhase();
875 876
    }
876 877

	
877 878
    /// \brief Runs the preflow algorithm to compute the minimum cut.
878 879
    ///
879 880
    /// Runs the preflow algorithm to compute the minimum cut.
880 881
    /// \note pf.runMinCut() is just a shortcut of the following code.
881 882
    /// \code
882 883
    ///   pf.init();
883 884
    ///   pf.startFirstPhase();
884 885
    /// \endcode
885 886
    void runMinCut() {
886 887
      init();
887 888
      startFirstPhase();
888 889
    }
889 890

	
890 891
    /// @}
891 892

	
892 893
    /// \name Query Functions
893 894
    /// The results of the preflow algorithm can be obtained using these
894 895
    /// functions.\n
895 896
    /// Either one of the \ref run() "run*()" functions or one of the
896 897
    /// \ref startFirstPhase() "start*()" functions should be called
897 898
    /// before using them.
898 899

	
899 900
    ///@{
900 901

	
901 902
    /// \brief Returns the value of the maximum flow.
902 903
    ///
903 904
    /// Returns the value of the maximum flow by returning the excess
904 905
    /// of the target node. This value equals to the value of
905 906
    /// the maximum flow already after the first phase of the algorithm.
906 907
    ///
907 908
    /// \pre Either \ref run() or \ref init() must be called before
908 909
    /// using this function.
909 910
    Value flowValue() const {
910 911
      return (*_excess)[_target];
911 912
    }
912 913

	
913 914
    /// \brief Returns the flow value on the given arc.
914 915
    ///
915 916
    /// Returns the flow value on the given arc. This method can
916 917
    /// be called after the second phase of the algorithm.
917 918
    ///
918 919
    /// \pre Either \ref run() or \ref init() must be called before
919 920
    /// using this function.
920 921
    Value flow(const Arc& arc) const {
921 922
      return (*_flow)[arc];
922 923
    }
923 924

	
924 925
    /// \brief Returns a const reference to the flow map.
925 926
    ///
926 927
    /// Returns a const reference to the arc map storing the found flow.
927 928
    /// This method can be called after the second phase of the algorithm.
928 929
    ///
929 930
    /// \pre Either \ref run() or \ref init() must be called before
930 931
    /// using this function.
931 932
    const FlowMap& flowMap() const {
932 933
      return *_flow;
933 934
    }
934 935

	
935 936
    /// \brief Returns \c true when the node is on the source side of the
936 937
    /// minimum cut.
937 938
    ///
938 939
    /// Returns true when the node is on the source side of the found
939 940
    /// minimum cut. This method can be called both after running \ref
940 941
    /// startFirstPhase() and \ref startSecondPhase().
941 942
    ///
942 943
    /// \pre Either \ref run() or \ref init() must be called before
943 944
    /// using this function.
944 945
    bool minCut(const Node& node) const {
945 946
      return ((*_level)[node] == _level->maxLevel()) == _phase;
946 947
    }
947 948

	
948 949
    /// \brief Gives back a minimum value cut.
949 950
    ///
950 951
    /// Sets \c cutMap to the characteristic vector of a minimum value
951 952
    /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
952 953
    /// node map with \c bool (or convertible) value type.
953 954
    ///
954 955
    /// This method can be called both after running \ref startFirstPhase()
955 956
    /// and \ref startSecondPhase(). The result after the second phase
956 957
    /// could be slightly different if inexact computation is used.
957 958
    ///
958 959
    /// \note This function calls \ref minCut() for each node, so it runs in
959 960
    /// O(n) time.
960 961
    ///
961 962
    /// \pre Either \ref run() or \ref init() must be called before
962 963
    /// using this function.
963 964
    template <typename CutMap>
964 965
    void minCutMap(CutMap& cutMap) const {
965 966
      for (NodeIt n(_graph); n != INVALID; ++n) {
966 967
        cutMap.set(n, minCut(n));
967 968
      }
968 969
    }
969 970

	
970 971
    /// @}
971 972
  };
972 973
}
973 974

	
974 975
#endif
0 comments (0 inline)