gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Add citations to the scaling MCF algorithms (#180, #184) and improve the doc of their group.
0 3 0
default
3 files changed with 12 insertions and 11 deletions:
↑ Collapse diff ↑
Show white space 1536 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 319
\e breadth-first \e search (BFS) and \e depth-first \e search (DFS)
320 320
\ref clrs01algorithms.
321 321
*/
322 322

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

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

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

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

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

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

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

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

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

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

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

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

	
395 395
/**
396 396
@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms
397 397
@ingroup algs
398 398

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

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

	
406 406
LEMON contains several algorithms for this problem.
407 407
 - \ref NetworkSimplex Primal Network Simplex algorithm with various
408 408
   pivot strategies \ref dantzig63linearprog, \ref kellyoneill91netsimplex.
409
 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
410
   cost scaling \ref goldberg90approximation, \ref goldberg97efficient,
409
 - \ref CostScaling Cost Scaling algorithm based on push/augment and
410
   relabel operations \ref goldberg90approximation, \ref goldberg97efficient,
411 411
   \ref bunnagel98efficient.
412
 - \ref CapacityScaling Successive Shortest %Path algorithm with optional
413
   capacity scaling \ref edmondskarp72theoretical.
414
 - \ref CancelAndTighten The Cancel and Tighten algorithm
415
   \ref goldberg89cyclecanceling.
416
 - \ref CycleCanceling Cycle-Canceling algorithms
417
   \ref klein67primal, \ref goldberg89cyclecanceling.
412
 - \ref CapacityScaling Capacity Scaling algorithm based on the successive
413
   shortest path method \ref edmondskarp72theoretical.
414
 - \ref CycleCanceling Cycle-Canceling algorithms, two of which are
415
   strongly polynomial \ref klein67primal, \ref goldberg89cyclecanceling.
418 416

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

	
425 423
/**
426 424
@defgroup min_cut Minimum Cut Algorithms
427 425
@ingroup algs
428 426

	
429 427
\brief Algorithms for finding minimum cut in graphs.
430 428

	
431 429
This group contains the algorithms for finding minimum cut in graphs.
432 430

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

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

	
442 440
LEMON contains several algorithms related to minimum cut problems:
443 441

	
444 442
- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut
445 443
  in directed graphs.
446 444
- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for
447 445
  calculating minimum cut in undirected graphs.
448 446
- \ref GomoryHu "Gomory-Hu tree computation" for calculating
449 447
  all-pairs minimum cut in undirected graphs.
450 448

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

	
455 453
/**
456 454
@defgroup min_mean_cycle Minimum Mean Cycle Algorithms
457 455
@ingroup algs
458 456
\brief Algorithms for finding minimum mean cycles.
459 457

	
460 458
This group contains the algorithms for finding minimum mean cycles
461 459
\ref clrs01algorithms, \ref amo93networkflows.
462 460

	
463 461
The \e minimum \e mean \e cycle \e problem is to find a directed cycle
464 462
of minimum mean length (cost) in a digraph.
465 463
The mean length of a cycle is the average length of its arcs, i.e. the
466 464
ratio between the total length of the cycle and the number of arcs on it.
467 465

	
468 466
This problem has an important connection to \e conservative \e length
469 467
\e functions, too. A length function on the arcs of a digraph is called
470 468
conservative if and only if there is no directed cycle of negative total
471 469
length. For an arbitrary length function, the negative of the minimum
472 470
cycle mean is the smallest \f$\epsilon\f$ value so that increasing the
473 471
arc lengths uniformly by \f$\epsilon\f$ results in a conservative length
474 472
function.
475 473

	
476 474
LEMON contains three algorithms for solving the minimum mean cycle problem:
477 475
- \ref Karp "Karp"'s original algorithm \ref amo93networkflows,
478 476
  \ref dasdan98minmeancycle.
479 477
- \ref HartmannOrlin "Hartmann-Orlin"'s algorithm, which is an improved
480 478
  version of Karp's algorithm \ref dasdan98minmeancycle.
481 479
- \ref Howard "Howard"'s policy iteration algorithm
482 480
  \ref dasdan98minmeancycle.
483 481

	
484 482
In practice, the Howard algorithm proved to be by far the most efficient
485 483
one, though the best known theoretical bound on its running time is
486 484
exponential.
487 485
Both Karp and HartmannOrlin algorithms run in time O(ne) and use space
488 486
O(n<sup>2</sup>+e), but the latter one is typically faster due to the
489 487
applied early termination scheme.
490 488
*/
491 489

	
492 490
/**
493 491
@defgroup matching Matching Algorithms
494 492
@ingroup algs
495 493
\brief Algorithms for finding matchings in graphs and bipartite graphs.
496 494

	
497 495
This group contains the algorithms for calculating
498 496
matchings in graphs and bipartite graphs. The general matching problem is
499 497
finding a subset of the edges for which each node has at most one incident
500 498
edge.
501 499

	
502 500
There are several different algorithms for calculate matchings in
503 501
graphs.  The matching problems in bipartite graphs are generally
504 502
easier than in general graphs. The goal of the matching optimization
505 503
can be finding maximum cardinality, maximum weight or minimum cost
506 504
matching. The search can be constrained to find perfect or
507 505
maximum cardinality matching.
508 506

	
509 507
The matching algorithms implemented in LEMON:
510 508
- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm
511 509
  for calculating maximum cardinality matching in bipartite graphs.
512 510
- \ref PrBipartiteMatching Push-relabel algorithm
513 511
  for calculating maximum cardinality matching in bipartite graphs.
514 512
- \ref MaxWeightedBipartiteMatching
515 513
  Successive shortest path algorithm for calculating maximum weighted
516 514
  matching and maximum weighted bipartite matching in bipartite graphs.
517 515
- \ref MinCostMaxBipartiteMatching
518 516
  Successive shortest path algorithm for calculating minimum cost maximum
519 517
  matching in bipartite graphs.
520 518
- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating
521 519
  maximum cardinality matching in general graphs.
522 520
- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating
523 521
  maximum weighted matching in general graphs.
524 522
- \ref MaxWeightedPerfectMatching
525 523
  Edmond's blossom shrinking algorithm for calculating maximum weighted
526 524
  perfect matching in general graphs.
527 525

	
528 526
\image html bipartite_matching.png
529 527
\image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth
530 528
*/
531 529

	
532 530
/**
533 531
@defgroup graph_properties Connectivity and Other Graph Properties
534 532
@ingroup algs
535 533
\brief Algorithms for discovering the graph properties
536 534

	
537 535
This group contains the algorithms for discovering the graph properties
538 536
like connectivity, bipartiteness, euler property, simplicity etc.
539 537

	
540 538
\image html connected_components.png
541 539
\image latex connected_components.eps "Connected components" width=\textwidth
542 540
*/
543 541

	
544 542
/**
545 543
@defgroup planar Planarity Embedding and Drawing
546 544
@ingroup algs
547 545
\brief Algorithms for planarity checking, embedding and drawing
548 546

	
549 547
This group contains the algorithms for planarity checking,
550 548
embedding and drawing.
551 549

	
552 550
\image html planar.png
553 551
\image latex planar.eps "Plane graph" width=\textwidth
554 552
*/
555 553

	
556 554
/**
557 555
@defgroup approx Approximation Algorithms
558 556
@ingroup algs
559 557
\brief Approximation algorithms.
560 558

	
561 559
This group contains the approximation and heuristic algorithms
562 560
implemented in LEMON.
563 561
*/
564 562

	
565 563
/**
566 564
@defgroup auxalg Auxiliary Algorithms
567 565
@ingroup algs
568 566
\brief Auxiliary algorithms implemented in LEMON.
569 567

	
570 568
This group contains some algorithms implemented in LEMON
571 569
in order to make it easier to implement complex algorithms.
572 570
*/
573 571

	
574 572
/**
575 573
@defgroup gen_opt_group General Optimization Tools
576 574
\brief This group contains some general optimization frameworks
577 575
implemented in LEMON.
578 576

	
579 577
This group contains some general optimization frameworks
580 578
implemented in LEMON.
581 579
*/
582 580

	
583 581
/**
584 582
@defgroup lp_group LP and MIP Solvers
585 583
@ingroup gen_opt_group
586 584
\brief LP and MIP solver interfaces for LEMON.
587 585

	
588 586
This group contains LP and MIP solver interfaces for LEMON.
589 587
Various LP solvers could be used in the same manner with this
590 588
high-level interface.
591 589

	
592 590
The currently supported solvers are \ref glpk, \ref clp, \ref cbc,
593 591
\ref cplex, \ref soplex.
594 592
*/
595 593

	
596 594
/**
597 595
@defgroup lp_utils Tools for Lp and Mip Solvers
598 596
@ingroup lp_group
599 597
\brief Helper tools to the Lp and Mip solvers.
600 598

	
601 599
This group adds some helper tools to general optimization framework
602 600
implemented in LEMON.
603 601
*/
604 602

	
605 603
/**
606 604
@defgroup metah Metaheuristics
607 605
@ingroup gen_opt_group
608 606
\brief Metaheuristics for LEMON library.
609 607

	
610 608
This group contains some metaheuristic optimization tools.
611 609
*/
612 610

	
613 611
/**
614 612
@defgroup utils Tools and Utilities
615 613
\brief Tools and utilities for programming in LEMON
616 614

	
617 615
Tools and utilities for programming in LEMON.
618 616
*/
619 617

	
620 618
/**
621 619
@defgroup gutils Basic Graph Utilities
622 620
@ingroup utils
623 621
\brief Simple basic graph utilities.
624 622

	
625 623
This group contains some simple basic graph utilities.
626 624
*/
627 625

	
628 626
/**
629 627
@defgroup misc Miscellaneous Tools
630 628
@ingroup utils
631 629
\brief Tools for development, debugging and testing.
632 630

	
633 631
This group contains several useful tools for development,
634 632
debugging and testing.
635 633
*/
636 634

	
637 635
/**
638 636
@defgroup timecount Time Measuring and Counting
639 637
@ingroup misc
640 638
\brief Simple tools for measuring the performance of algorithms.
641 639

	
642 640
This group contains simple tools for measuring the performance
643 641
of algorithms.
644 642
*/
645 643

	
646 644
/**
647 645
@defgroup exceptions Exceptions
648 646
@ingroup utils
649 647
\brief Exceptions defined in LEMON.
650 648

	
651 649
This group contains the exceptions defined in LEMON.
652 650
*/
653 651

	
654 652
/**
655 653
@defgroup io_group Input-Output
656 654
\brief Graph Input-Output methods
657 655

	
658 656
This group contains the tools for importing and exporting graphs
659 657
and graph related data. Now it supports the \ref lgf-format
660 658
"LEMON Graph Format", the \c DIMACS format and the encapsulated
661 659
postscript (EPS) format.
662 660
*/
663 661

	
664 662
/**
665 663
@defgroup lemon_io LEMON Graph Format
666 664
@ingroup io_group
667 665
\brief Reading and writing LEMON Graph Format.
668 666

	
669 667
This group contains methods for reading and writing
670 668
\ref lgf-format "LEMON Graph Format".
671 669
*/
672 670

	
673 671
/**
674 672
@defgroup eps_io Postscript Exporting
675 673
@ingroup io_group
676 674
\brief General \c EPS drawer and graph exporter
677 675

	
678 676
This group contains general \c EPS drawing methods and special
679 677
graph exporting tools.
680 678
*/
681 679

	
682 680
/**
683 681
@defgroup dimacs_group DIMACS Format
684 682
@ingroup io_group
685 683
\brief Read and write files in DIMACS format
686 684

	
687 685
Tools to read a digraph from or write it to a file in DIMACS format data.
688 686
*/
689 687

	
690 688
/**
691 689
@defgroup nauty_group NAUTY Format
692 690
@ingroup io_group
693 691
\brief Read \e Nauty format
694 692

	
695 693
Tool to read graphs from \e Nauty format data.
696 694
*/
697 695

	
698 696
/**
699 697
@defgroup concept Concepts
700 698
\brief Skeleton classes and concept checking classes
701 699

	
702 700
This group contains the data/algorithm skeletons and concept checking
703 701
classes implemented in LEMON.
704 702

	
705 703
The purpose of the classes in this group is fourfold.
706 704

	
707 705
- These classes contain the documentations of the %concepts. In order
708 706
  to avoid document multiplications, an implementation of a concept
709 707
  simply refers to the corresponding concept class.
710 708

	
711 709
- These classes declare every functions, <tt>typedef</tt>s etc. an
712 710
  implementation of the %concepts should provide, however completely
713 711
  without implementations and real data structures behind the
714 712
  interface. On the other hand they should provide nothing else. All
715 713
  the algorithms working on a data structure meeting a certain concept
716 714
  should compile with these classes. (Though it will not run properly,
717 715
  of course.) In this way it is easily to check if an algorithm
718 716
  doesn't use any extra feature of a certain implementation.
719 717

	
720 718
- The concept descriptor classes also provide a <em>checker class</em>
721 719
  that makes it possible to check whether a certain implementation of a
722 720
  concept indeed provides all the required features.
723 721

	
724 722
- Finally, They can serve as a skeleton of a new implementation of a concept.
725 723
*/
726 724

	
727 725
/**
728 726
@defgroup graph_concepts Graph Structure Concepts
729 727
@ingroup concept
730 728
\brief Skeleton and concept checking classes for graph structures
731 729

	
732 730
This group contains the skeletons and concept checking classes of
733 731
graph structures.
734 732
*/
735 733

	
736 734
/**
737 735
@defgroup map_concepts Map Concepts
738 736
@ingroup concept
739 737
\brief Skeleton and concept checking classes for maps
740 738

	
741 739
This group contains the skeletons and concept checking classes of maps.
742 740
*/
743 741

	
744 742
/**
745 743
@defgroup tools Standalone Utility Applications
746 744

	
747 745
Some utility applications are listed here.
748 746

	
749 747
The standard compilation procedure (<tt>./configure;make</tt>) will compile
750 748
them, as well.
751 749
*/
752 750

	
753 751
/**
754 752
\anchor demoprograms
755 753

	
756 754
@defgroup demos Demo Programs
757 755

	
758 756
Some demo programs are listed here. Their full source codes can be found in
759 757
the \c demo subdirectory of the source tree.
760 758

	
761 759
In order to compile them, use the <tt>make demo</tt> or the
762 760
<tt>make check</tt> commands.
763 761
*/
764 762

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

	
19 19
#ifndef LEMON_CAPACITY_SCALING_H
20 20
#define LEMON_CAPACITY_SCALING_H
21 21

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

	
27 27
#include <vector>
28 28
#include <limits>
29 29
#include <lemon/core.h>
30 30
#include <lemon/bin_heap.h>
31 31

	
32 32
namespace lemon {
33 33

	
34 34
  /// \brief Default traits class of CapacityScaling algorithm.
35 35
  ///
36 36
  /// Default traits class of CapacityScaling algorithm.
37 37
  /// \tparam GR Digraph type.
38 38
  /// \tparam V The number type used for flow amounts, capacity bounds
39 39
  /// and supply values. By default it is \c int.
40 40
  /// \tparam C The number type used for costs and potentials.
41 41
  /// By default it is the same as \c V.
42 42
  template <typename GR, typename V = int, typename C = V>
43 43
  struct CapacityScalingDefaultTraits
44 44
  {
45 45
    /// The type of the digraph
46 46
    typedef GR Digraph;
47 47
    /// The type of the flow amounts, capacity bounds and supply values
48 48
    typedef V Value;
49 49
    /// The type of the arc costs
50 50
    typedef C Cost;
51 51

	
52 52
    /// \brief The type of the heap used for internal Dijkstra computations.
53 53
    ///
54 54
    /// The type of the heap used for internal Dijkstra computations.
55 55
    /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
56 56
    /// its priority type must be \c Cost and its cross reference type
57 57
    /// must be \ref RangeMap "RangeMap<int>".
58 58
    typedef BinHeap<Cost, RangeMap<int> > Heap;
59 59
  };
60 60

	
61 61
  /// \addtogroup min_cost_flow_algs
62 62
  /// @{
63 63

	
64 64
  /// \brief Implementation of the Capacity Scaling algorithm for
65 65
  /// finding a \ref min_cost_flow "minimum cost flow".
66 66
  ///
67 67
  /// \ref CapacityScaling implements the capacity scaling version
68 68
  /// of the successive shortest path algorithm for finding a
69
  /// \ref min_cost_flow "minimum cost flow". It is an efficient dual
69
  /// \ref min_cost_flow "minimum cost flow" \ref amo93networkflows,
70
  /// \ref edmondskarp72theoretical. It is an efficient dual
70 71
  /// solution method.
71 72
  ///
72 73
  /// Most of the parameters of the problem (except for the digraph)
73 74
  /// can be given using separate functions, and the algorithm can be
74 75
  /// executed using the \ref run() function. If some parameters are not
75 76
  /// specified, then default values will be used.
76 77
  ///
77 78
  /// \tparam GR The digraph type the algorithm runs on.
78 79
  /// \tparam V The number type used for flow amounts, capacity bounds
79 80
  /// and supply values in the algorithm. By default it is \c int.
80 81
  /// \tparam C The number type used for costs and potentials in the
81 82
  /// algorithm. By default it is the same as \c V.
82 83
  ///
83 84
  /// \warning Both number types must be signed and all input data must
84 85
  /// be integer.
85 86
  /// \warning This algorithm does not support negative costs for such
86 87
  /// arcs that have infinite upper bound.
87 88
#ifdef DOXYGEN
88 89
  template <typename GR, typename V, typename C, typename TR>
89 90
#else
90 91
  template < typename GR, typename V = int, typename C = V,
91 92
             typename TR = CapacityScalingDefaultTraits<GR, V, C> >
92 93
#endif
93 94
  class CapacityScaling
94 95
  {
95 96
  public:
96 97

	
97 98
    /// The type of the digraph
98 99
    typedef typename TR::Digraph Digraph;
99 100
    /// The type of the flow amounts, capacity bounds and supply values
100 101
    typedef typename TR::Value Value;
101 102
    /// The type of the arc costs
102 103
    typedef typename TR::Cost Cost;
103 104

	
104 105
    /// The type of the heap used for internal Dijkstra computations
105 106
    typedef typename TR::Heap Heap;
106 107

	
107 108
    /// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm
108 109
    typedef TR Traits;
109 110

	
110 111
  public:
111 112

	
112 113
    /// \brief Problem type constants for the \c run() function.
113 114
    ///
114 115
    /// Enum type containing the problem type constants that can be
115 116
    /// returned by the \ref run() function of the algorithm.
116 117
    enum ProblemType {
117 118
      /// The problem has no feasible solution (flow).
118 119
      INFEASIBLE,
119 120
      /// The problem has optimal solution (i.e. it is feasible and
120 121
      /// bounded), and the algorithm has found optimal flow and node
121 122
      /// potentials (primal and dual solutions).
122 123
      OPTIMAL,
123 124
      /// The digraph contains an arc of negative cost and infinite
124 125
      /// upper bound. It means that the objective function is unbounded
125 126
      /// on that arc, however, note that it could actually be bounded
126 127
      /// over the feasible flows, but this algroithm cannot handle
127 128
      /// these cases.
128 129
      UNBOUNDED
129 130
    };
130 131
  
131 132
  private:
132 133

	
133 134
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
134 135

	
135 136
    typedef std::vector<int> IntVector;
136 137
    typedef std::vector<char> BoolVector;
137 138
    typedef std::vector<Value> ValueVector;
138 139
    typedef std::vector<Cost> CostVector;
139 140

	
140 141
  private:
141 142

	
142 143
    // Data related to the underlying digraph
143 144
    const GR &_graph;
144 145
    int _node_num;
145 146
    int _arc_num;
146 147
    int _res_arc_num;
147 148
    int _root;
148 149

	
149 150
    // Parameters of the problem
150 151
    bool _have_lower;
151 152
    Value _sum_supply;
152 153

	
153 154
    // Data structures for storing the digraph
154 155
    IntNodeMap _node_id;
155 156
    IntArcMap _arc_idf;
156 157
    IntArcMap _arc_idb;
157 158
    IntVector _first_out;
158 159
    BoolVector _forward;
159 160
    IntVector _source;
160 161
    IntVector _target;
161 162
    IntVector _reverse;
162 163

	
163 164
    // Node and arc data
164 165
    ValueVector _lower;
165 166
    ValueVector _upper;
166 167
    CostVector _cost;
167 168
    ValueVector _supply;
168 169

	
169 170
    ValueVector _res_cap;
170 171
    CostVector _pi;
171 172
    ValueVector _excess;
172 173
    IntVector _excess_nodes;
173 174
    IntVector _deficit_nodes;
174 175

	
175 176
    Value _delta;
176 177
    int _factor;
177 178
    IntVector _pred;
178 179

	
179 180
  public:
180 181
  
181 182
    /// \brief Constant for infinite upper bounds (capacities).
182 183
    ///
183 184
    /// Constant for infinite upper bounds (capacities).
184 185
    /// It is \c std::numeric_limits<Value>::infinity() if available,
185 186
    /// \c std::numeric_limits<Value>::max() otherwise.
186 187
    const Value INF;
187 188

	
188 189
  private:
189 190

	
190 191
    // Special implementation of the Dijkstra algorithm for finding
191 192
    // shortest paths in the residual network of the digraph with
192 193
    // respect to the reduced arc costs and modifying the node
193 194
    // potentials according to the found distance labels.
194 195
    class ResidualDijkstra
195 196
    {
196 197
    private:
197 198

	
198 199
      int _node_num;
199 200
      bool _geq;
200 201
      const IntVector &_first_out;
201 202
      const IntVector &_target;
202 203
      const CostVector &_cost;
203 204
      const ValueVector &_res_cap;
204 205
      const ValueVector &_excess;
205 206
      CostVector &_pi;
206 207
      IntVector &_pred;
207 208
      
208 209
      IntVector _proc_nodes;
209 210
      CostVector _dist;
210 211
      
211 212
    public:
212 213

	
213 214
      ResidualDijkstra(CapacityScaling& cs) :
214 215
        _node_num(cs._node_num), _geq(cs._sum_supply < 0),
215 216
        _first_out(cs._first_out), _target(cs._target), _cost(cs._cost),
216 217
        _res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi),
217 218
        _pred(cs._pred), _dist(cs._node_num)
218 219
      {}
219 220

	
220 221
      int run(int s, Value delta = 1) {
221 222
        RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP);
222 223
        Heap heap(heap_cross_ref);
223 224
        heap.push(s, 0);
224 225
        _pred[s] = -1;
225 226
        _proc_nodes.clear();
226 227

	
227 228
        // Process nodes
228 229
        while (!heap.empty() && _excess[heap.top()] > -delta) {
229 230
          int u = heap.top(), v;
230 231
          Cost d = heap.prio() + _pi[u], dn;
231 232
          _dist[u] = heap.prio();
232 233
          _proc_nodes.push_back(u);
233 234
          heap.pop();
234 235

	
235 236
          // Traverse outgoing residual arcs
236 237
          int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1;
237 238
          for (int a = _first_out[u]; a != last_out; ++a) {
238 239
            if (_res_cap[a] < delta) continue;
239 240
            v = _target[a];
240 241
            switch (heap.state(v)) {
241 242
              case Heap::PRE_HEAP:
242 243
                heap.push(v, d + _cost[a] - _pi[v]);
243 244
                _pred[v] = a;
244 245
                break;
245 246
              case Heap::IN_HEAP:
246 247
                dn = d + _cost[a] - _pi[v];
247 248
                if (dn < heap[v]) {
248 249
                  heap.decrease(v, dn);
249 250
                  _pred[v] = a;
250 251
                }
251 252
                break;
252 253
              case Heap::POST_HEAP:
253 254
                break;
254 255
            }
255 256
          }
256 257
        }
257 258
        if (heap.empty()) return -1;
258 259

	
259 260
        // Update potentials of processed nodes
260 261
        int t = heap.top();
261 262
        Cost dt = heap.prio();
262 263
        for (int i = 0; i < int(_proc_nodes.size()); ++i) {
263 264
          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;
264 265
        }
265 266

	
266 267
        return t;
267 268
      }
268 269

	
269 270
    }; //class ResidualDijkstra
270 271

	
271 272
  public:
272 273

	
273 274
    /// \name Named Template Parameters
274 275
    /// @{
275 276

	
276 277
    template <typename T>
277 278
    struct SetHeapTraits : public Traits {
278 279
      typedef T Heap;
279 280
    };
280 281

	
281 282
    /// \brief \ref named-templ-param "Named parameter" for setting
282 283
    /// \c Heap type.
283 284
    ///
284 285
    /// \ref named-templ-param "Named parameter" for setting \c Heap
285 286
    /// type, which is used for internal Dijkstra computations.
286 287
    /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
287 288
    /// its priority type must be \c Cost and its cross reference type
288 289
    /// must be \ref RangeMap "RangeMap<int>".
289 290
    template <typename T>
290 291
    struct SetHeap
291 292
      : public CapacityScaling<GR, V, C, SetHeapTraits<T> > {
292 293
      typedef  CapacityScaling<GR, V, C, SetHeapTraits<T> > Create;
293 294
    };
294 295

	
295 296
    /// @}
296 297

	
297 298
  public:
298 299

	
299 300
    /// \brief Constructor.
300 301
    ///
301 302
    /// The constructor of the class.
302 303
    ///
303 304
    /// \param graph The digraph the algorithm runs on.
304 305
    CapacityScaling(const GR& graph) :
305 306
      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
306 307
      INF(std::numeric_limits<Value>::has_infinity ?
307 308
          std::numeric_limits<Value>::infinity() :
308 309
          std::numeric_limits<Value>::max())
309 310
    {
310 311
      // Check the number types
311 312
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
312 313
        "The flow type of CapacityScaling must be signed");
313 314
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
314 315
        "The cost type of CapacityScaling must be signed");
315 316

	
316 317
      // Resize vectors
317 318
      _node_num = countNodes(_graph);
318 319
      _arc_num = countArcs(_graph);
319 320
      _res_arc_num = 2 * (_arc_num + _node_num);
320 321
      _root = _node_num;
321 322
      ++_node_num;
322 323

	
323 324
      _first_out.resize(_node_num + 1);
324 325
      _forward.resize(_res_arc_num);
325 326
      _source.resize(_res_arc_num);
326 327
      _target.resize(_res_arc_num);
327 328
      _reverse.resize(_res_arc_num);
328 329

	
329 330
      _lower.resize(_res_arc_num);
330 331
      _upper.resize(_res_arc_num);
331 332
      _cost.resize(_res_arc_num);
332 333
      _supply.resize(_node_num);
333 334
      
334 335
      _res_cap.resize(_res_arc_num);
335 336
      _pi.resize(_node_num);
336 337
      _excess.resize(_node_num);
337 338
      _pred.resize(_node_num);
338 339

	
339 340
      // Copy the graph
340 341
      int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
341 342
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
342 343
        _node_id[n] = i;
343 344
      }
344 345
      i = 0;
345 346
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
346 347
        _first_out[i] = j;
347 348
        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
348 349
          _arc_idf[a] = j;
349 350
          _forward[j] = true;
350 351
          _source[j] = i;
351 352
          _target[j] = _node_id[_graph.runningNode(a)];
352 353
        }
353 354
        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
354 355
          _arc_idb[a] = j;
355 356
          _forward[j] = false;
356 357
          _source[j] = i;
357 358
          _target[j] = _node_id[_graph.runningNode(a)];
358 359
        }
359 360
        _forward[j] = false;
360 361
        _source[j] = i;
361 362
        _target[j] = _root;
362 363
        _reverse[j] = k;
363 364
        _forward[k] = true;
364 365
        _source[k] = _root;
365 366
        _target[k] = i;
366 367
        _reverse[k] = j;
367 368
        ++j; ++k;
368 369
      }
369 370
      _first_out[i] = j;
370 371
      _first_out[_node_num] = k;
371 372
      for (ArcIt a(_graph); a != INVALID; ++a) {
372 373
        int fi = _arc_idf[a];
373 374
        int bi = _arc_idb[a];
374 375
        _reverse[fi] = bi;
375 376
        _reverse[bi] = fi;
376 377
      }
377 378
      
378 379
      // Reset parameters
379 380
      reset();
380 381
    }
381 382

	
382 383
    /// \name Parameters
383 384
    /// The parameters of the algorithm can be specified using these
384 385
    /// functions.
385 386

	
386 387
    /// @{
387 388

	
388 389
    /// \brief Set the lower bounds on the arcs.
389 390
    ///
390 391
    /// This function sets the lower bounds on the arcs.
391 392
    /// If it is not used before calling \ref run(), the lower bounds
392 393
    /// will be set to zero on all arcs.
393 394
    ///
394 395
    /// \param map An arc map storing the lower bounds.
395 396
    /// Its \c Value type must be convertible to the \c Value type
396 397
    /// of the algorithm.
397 398
    ///
398 399
    /// \return <tt>(*this)</tt>
399 400
    template <typename LowerMap>
400 401
    CapacityScaling& lowerMap(const LowerMap& map) {
401 402
      _have_lower = true;
402 403
      for (ArcIt a(_graph); a != INVALID; ++a) {
403 404
        _lower[_arc_idf[a]] = map[a];
404 405
        _lower[_arc_idb[a]] = map[a];
405 406
      }
406 407
      return *this;
407 408
    }
408 409

	
409 410
    /// \brief Set the upper bounds (capacities) on the arcs.
410 411
    ///
411 412
    /// This function sets the upper bounds (capacities) on the arcs.
412 413
    /// If it is not used before calling \ref run(), the upper bounds
413 414
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
414 415
    /// unbounded from above).
415 416
    ///
416 417
    /// \param map An arc map storing the upper bounds.
417 418
    /// Its \c Value type must be convertible to the \c Value type
418 419
    /// of the algorithm.
419 420
    ///
420 421
    /// \return <tt>(*this)</tt>
421 422
    template<typename UpperMap>
422 423
    CapacityScaling& upperMap(const UpperMap& map) {
423 424
      for (ArcIt a(_graph); a != INVALID; ++a) {
424 425
        _upper[_arc_idf[a]] = map[a];
425 426
      }
426 427
      return *this;
427 428
    }
428 429

	
429 430
    /// \brief Set the costs of the arcs.
430 431
    ///
431 432
    /// This function sets the costs of the arcs.
432 433
    /// If it is not used before calling \ref run(), the costs
433 434
    /// will be set to \c 1 on all arcs.
434 435
    ///
435 436
    /// \param map An arc map storing the costs.
436 437
    /// Its \c Value type must be convertible to the \c Cost type
437 438
    /// of the algorithm.
438 439
    ///
439 440
    /// \return <tt>(*this)</tt>
440 441
    template<typename CostMap>
441 442
    CapacityScaling& costMap(const CostMap& map) {
442 443
      for (ArcIt a(_graph); a != INVALID; ++a) {
443 444
        _cost[_arc_idf[a]] =  map[a];
444 445
        _cost[_arc_idb[a]] = -map[a];
445 446
      }
446 447
      return *this;
447 448
    }
448 449

	
449 450
    /// \brief Set the supply values of the nodes.
450 451
    ///
451 452
    /// This function sets the supply values of the nodes.
452 453
    /// If neither this function nor \ref stSupply() is used before
453 454
    /// calling \ref run(), the supply of each node will be set to zero.
454 455
    ///
455 456
    /// \param map A node map storing the supply values.
456 457
    /// Its \c Value type must be convertible to the \c Value type
457 458
    /// of the algorithm.
458 459
    ///
459 460
    /// \return <tt>(*this)</tt>
460 461
    template<typename SupplyMap>
461 462
    CapacityScaling& supplyMap(const SupplyMap& map) {
462 463
      for (NodeIt n(_graph); n != INVALID; ++n) {
463 464
        _supply[_node_id[n]] = map[n];
464 465
      }
465 466
      return *this;
466 467
    }
467 468

	
468 469
    /// \brief Set single source and target nodes and a supply value.
469 470
    ///
470 471
    /// This function sets a single source node and a single target node
471 472
    /// and the required flow value.
472 473
    /// If neither this function nor \ref supplyMap() is used before
473 474
    /// calling \ref run(), the supply of each node will be set to zero.
474 475
    ///
475 476
    /// Using this function has the same effect as using \ref supplyMap()
476 477
    /// with such a map in which \c k is assigned to \c s, \c -k is
477 478
    /// assigned to \c t and all other nodes have zero supply value.
478 479
    ///
479 480
    /// \param s The source node.
480 481
    /// \param t The target node.
481 482
    /// \param k The required amount of flow from node \c s to node \c t
482 483
    /// (i.e. the supply of \c s and the demand of \c t).
483 484
    ///
484 485
    /// \return <tt>(*this)</tt>
485 486
    CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
486 487
      for (int i = 0; i != _node_num; ++i) {
487 488
        _supply[i] = 0;
488 489
      }
489 490
      _supply[_node_id[s]] =  k;
490 491
      _supply[_node_id[t]] = -k;
491 492
      return *this;
492 493
    }
493 494
    
494 495
    /// @}
495 496

	
496 497
    /// \name Execution control
497 498
    /// The algorithm can be executed using \ref run().
498 499

	
499 500
    /// @{
500 501

	
501 502
    /// \brief Run the algorithm.
502 503
    ///
503 504
    /// This function runs the algorithm.
504 505
    /// The paramters can be specified using functions \ref lowerMap(),
505 506
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
506 507
    /// For example,
507 508
    /// \code
508 509
    ///   CapacityScaling<ListDigraph> cs(graph);
509 510
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
510 511
    ///     .supplyMap(sup).run();
511 512
    /// \endcode
512 513
    ///
513 514
    /// This function can be called more than once. All the parameters
514 515
    /// that have been given are kept for the next call, unless
515 516
    /// \ref reset() is called, thus only the modified parameters
516 517
    /// have to be set again. See \ref reset() for examples.
517 518
    /// However, the underlying digraph must not be modified after this
518 519
    /// class have been constructed, since it copies and extends the graph.
519 520
    ///
520 521
    /// \param factor The capacity scaling factor. It must be larger than
521 522
    /// one to use scaling. If it is less or equal to one, then scaling
522 523
    /// will be disabled.
523 524
    ///
524 525
    /// \return \c INFEASIBLE if no feasible flow exists,
525 526
    /// \n \c OPTIMAL if the problem has optimal solution
526 527
    /// (i.e. it is feasible and bounded), and the algorithm has found
527 528
    /// optimal flow and node potentials (primal and dual solutions),
528 529
    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
529 530
    /// and infinite upper bound. It means that the objective function
530 531
    /// is unbounded on that arc, however, note that it could actually be
531 532
    /// bounded over the feasible flows, but this algroithm cannot handle
532 533
    /// these cases.
533 534
    ///
534 535
    /// \see ProblemType
535 536
    ProblemType run(int factor = 4) {
536 537
      _factor = factor;
537 538
      ProblemType pt = init();
538 539
      if (pt != OPTIMAL) return pt;
539 540
      return start();
540 541
    }
541 542

	
542 543
    /// \brief Reset all the parameters that have been given before.
543 544
    ///
544 545
    /// This function resets all the paramaters that have been given
545 546
    /// before using functions \ref lowerMap(), \ref upperMap(),
546 547
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
547 548
    ///
548 549
    /// It is useful for multiple run() calls. If this function is not
549 550
    /// used, all the parameters given before are kept for the next
550 551
    /// \ref run() call.
551 552
    /// However, the underlying digraph must not be modified after this
552 553
    /// class have been constructed, since it copies and extends the graph.
553 554
    ///
554 555
    /// For example,
555 556
    /// \code
556 557
    ///   CapacityScaling<ListDigraph> cs(graph);
557 558
    ///
558 559
    ///   // First run
559 560
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
560 561
    ///     .supplyMap(sup).run();
561 562
    ///
562 563
    ///   // Run again with modified cost map (reset() is not called,
563 564
    ///   // so only the cost map have to be set again)
564 565
    ///   cost[e] += 100;
565 566
    ///   cs.costMap(cost).run();
566 567
    ///
567 568
    ///   // Run again from scratch using reset()
568 569
    ///   // (the lower bounds will be set to zero on all arcs)
569 570
    ///   cs.reset();
570 571
    ///   cs.upperMap(capacity).costMap(cost)
571 572
    ///     .supplyMap(sup).run();
572 573
    /// \endcode
573 574
    ///
574 575
    /// \return <tt>(*this)</tt>
575 576
    CapacityScaling& reset() {
576 577
      for (int i = 0; i != _node_num; ++i) {
577 578
        _supply[i] = 0;
578 579
      }
579 580
      for (int j = 0; j != _res_arc_num; ++j) {
580 581
        _lower[j] = 0;
581 582
        _upper[j] = INF;
582 583
        _cost[j] = _forward[j] ? 1 : -1;
583 584
      }
584 585
      _have_lower = false;
585 586
      return *this;
586 587
    }
587 588

	
588 589
    /// @}
589 590

	
590 591
    /// \name Query Functions
591 592
    /// The results of the algorithm can be obtained using these
592 593
    /// functions.\n
593 594
    /// The \ref run() function must be called before using them.
594 595

	
595 596
    /// @{
596 597

	
597 598
    /// \brief Return the total cost of the found flow.
598 599
    ///
599 600
    /// This function returns the total cost of the found flow.
600 601
    /// Its complexity is O(e).
601 602
    ///
602 603
    /// \note The return type of the function can be specified as a
603 604
    /// template parameter. For example,
604 605
    /// \code
605 606
    ///   cs.totalCost<double>();
606 607
    /// \endcode
607 608
    /// It is useful if the total cost cannot be stored in the \c Cost
608 609
    /// type of the algorithm, which is the default return type of the
609 610
    /// function.
610 611
    ///
611 612
    /// \pre \ref run() must be called before using this function.
612 613
    template <typename Number>
613 614
    Number totalCost() const {
614 615
      Number c = 0;
615 616
      for (ArcIt a(_graph); a != INVALID; ++a) {
616 617
        int i = _arc_idb[a];
617 618
        c += static_cast<Number>(_res_cap[i]) *
618 619
             (-static_cast<Number>(_cost[i]));
619 620
      }
620 621
      return c;
621 622
    }
622 623

	
623 624
#ifndef DOXYGEN
624 625
    Cost totalCost() const {
625 626
      return totalCost<Cost>();
626 627
    }
627 628
#endif
628 629

	
629 630
    /// \brief Return the flow on the given arc.
630 631
    ///
631 632
    /// This function returns the flow on the given arc.
632 633
    ///
633 634
    /// \pre \ref run() must be called before using this function.
634 635
    Value flow(const Arc& a) const {
635 636
      return _res_cap[_arc_idb[a]];
636 637
    }
637 638

	
638 639
    /// \brief Return the flow map (the primal solution).
639 640
    ///
640 641
    /// This function copies the flow value on each arc into the given
641 642
    /// map. The \c Value type of the algorithm must be convertible to
642 643
    /// the \c Value type of the map.
643 644
    ///
644 645
    /// \pre \ref run() must be called before using this function.
645 646
    template <typename FlowMap>
646 647
    void flowMap(FlowMap &map) const {
647 648
      for (ArcIt a(_graph); a != INVALID; ++a) {
648 649
        map.set(a, _res_cap[_arc_idb[a]]);
649 650
      }
650 651
    }
651 652

	
652 653
    /// \brief Return the potential (dual value) of the given node.
653 654
    ///
654 655
    /// This function returns the potential (dual value) of the
655 656
    /// given node.
656 657
    ///
657 658
    /// \pre \ref run() must be called before using this function.
658 659
    Cost potential(const Node& n) const {
659 660
      return _pi[_node_id[n]];
660 661
    }
661 662

	
662 663
    /// \brief Return the potential map (the dual solution).
663 664
    ///
664 665
    /// This function copies the potential (dual value) of each node
665 666
    /// into the given map.
666 667
    /// The \c Cost type of the algorithm must be convertible to the
667 668
    /// \c Value type of the map.
668 669
    ///
669 670
    /// \pre \ref run() must be called before using this function.
670 671
    template <typename PotentialMap>
671 672
    void potentialMap(PotentialMap &map) const {
672 673
      for (NodeIt n(_graph); n != INVALID; ++n) {
673 674
        map.set(n, _pi[_node_id[n]]);
674 675
      }
675 676
    }
676 677

	
677 678
    /// @}
678 679

	
679 680
  private:
680 681

	
681 682
    // Initialize the algorithm
682 683
    ProblemType init() {
683 684
      if (_node_num == 0) return INFEASIBLE;
684 685

	
685 686
      // Check the sum of supply values
686 687
      _sum_supply = 0;
687 688
      for (int i = 0; i != _root; ++i) {
688 689
        _sum_supply += _supply[i];
689 690
      }
690 691
      if (_sum_supply > 0) return INFEASIBLE;
691 692
      
692 693
      // Initialize vectors
693 694
      for (int i = 0; i != _root; ++i) {
694 695
        _pi[i] = 0;
695 696
        _excess[i] = _supply[i];
696 697
      }
697 698

	
698 699
      // Remove non-zero lower bounds
699 700
      const Value MAX = std::numeric_limits<Value>::max();
700 701
      int last_out;
701 702
      if (_have_lower) {
702 703
        for (int i = 0; i != _root; ++i) {
703 704
          last_out = _first_out[i+1];
704 705
          for (int j = _first_out[i]; j != last_out; ++j) {
705 706
            if (_forward[j]) {
706 707
              Value c = _lower[j];
707 708
              if (c >= 0) {
708 709
                _res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF;
709 710
              } else {
710 711
                _res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF;
711 712
              }
712 713
              _excess[i] -= c;
713 714
              _excess[_target[j]] += c;
714 715
            } else {
715 716
              _res_cap[j] = 0;
716 717
            }
717 718
          }
718 719
        }
719 720
      } else {
720 721
        for (int j = 0; j != _res_arc_num; ++j) {
721 722
          _res_cap[j] = _forward[j] ? _upper[j] : 0;
722 723
        }
723 724
      }
724 725

	
725 726
      // Handle negative costs
726 727
      for (int i = 0; i != _root; ++i) {
727 728
        last_out = _first_out[i+1] - 1;
728 729
        for (int j = _first_out[i]; j != last_out; ++j) {
729 730
          Value rc = _res_cap[j];
730 731
          if (_cost[j] < 0 && rc > 0) {
731 732
            if (rc >= MAX) return UNBOUNDED;
732 733
            _excess[i] -= rc;
733 734
            _excess[_target[j]] += rc;
734 735
            _res_cap[j] = 0;
735 736
            _res_cap[_reverse[j]] += rc;
736 737
          }
737 738
        }
738 739
      }
739 740
      
740 741
      // Handle GEQ supply type
741 742
      if (_sum_supply < 0) {
742 743
        _pi[_root] = 0;
743 744
        _excess[_root] = -_sum_supply;
744 745
        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
745 746
          int ra = _reverse[a];
746 747
          _res_cap[a] = -_sum_supply + 1;
747 748
          _res_cap[ra] = 0;
748 749
          _cost[a] = 0;
749 750
          _cost[ra] = 0;
750 751
        }
751 752
      } else {
752 753
        _pi[_root] = 0;
753 754
        _excess[_root] = 0;
754 755
        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
755 756
          int ra = _reverse[a];
756 757
          _res_cap[a] = 1;
757 758
          _res_cap[ra] = 0;
758 759
          _cost[a] = 0;
759 760
          _cost[ra] = 0;
760 761
        }
761 762
      }
762 763

	
763 764
      // Initialize delta value
764 765
      if (_factor > 1) {
765 766
        // With scaling
766 767
        Value max_sup = 0, max_dem = 0;
767 768
        for (int i = 0; i != _node_num; ++i) {
768 769
          Value ex = _excess[i];
769 770
          if ( ex > max_sup) max_sup =  ex;
770 771
          if (-ex > max_dem) max_dem = -ex;
771 772
        }
772 773
        Value max_cap = 0;
773 774
        for (int j = 0; j != _res_arc_num; ++j) {
774 775
          if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
775 776
        }
776 777
        max_sup = std::min(std::min(max_sup, max_dem), max_cap);
777 778
        for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;
778 779
      } else {
779 780
        // Without scaling
780 781
        _delta = 1;
781 782
      }
782 783

	
783 784
      return OPTIMAL;
784 785
    }
785 786

	
786 787
    ProblemType start() {
787 788
      // Execute the algorithm
788 789
      ProblemType pt;
789 790
      if (_delta > 1)
790 791
        pt = startWithScaling();
791 792
      else
792 793
        pt = startWithoutScaling();
793 794

	
794 795
      // Handle non-zero lower bounds
795 796
      if (_have_lower) {
796 797
        int limit = _first_out[_root];
797 798
        for (int j = 0; j != limit; ++j) {
798 799
          if (!_forward[j]) _res_cap[j] += _lower[j];
799 800
        }
800 801
      }
801 802

	
802 803
      // Shift potentials if necessary
803 804
      Cost pr = _pi[_root];
804 805
      if (_sum_supply < 0 || pr > 0) {
805 806
        for (int i = 0; i != _node_num; ++i) {
806 807
          _pi[i] -= pr;
807 808
        }        
808 809
      }
809 810
      
810 811
      return pt;
811 812
    }
812 813

	
813 814
    // Execute the capacity scaling algorithm
814 815
    ProblemType startWithScaling() {
815 816
      // Perform capacity scaling phases
816 817
      int s, t;
817 818
      ResidualDijkstra _dijkstra(*this);
818 819
      while (true) {
819 820
        // Saturate all arcs not satisfying the optimality condition
820 821
        int last_out;
821 822
        for (int u = 0; u != _node_num; ++u) {
822 823
          last_out = _sum_supply < 0 ?
823 824
            _first_out[u+1] : _first_out[u+1] - 1;
824 825
          for (int a = _first_out[u]; a != last_out; ++a) {
825 826
            int v = _target[a];
826 827
            Cost c = _cost[a] + _pi[u] - _pi[v];
827 828
            Value rc = _res_cap[a];
828 829
            if (c < 0 && rc >= _delta) {
829 830
              _excess[u] -= rc;
830 831
              _excess[v] += rc;
831 832
              _res_cap[a] = 0;
832 833
              _res_cap[_reverse[a]] += rc;
833 834
            }
834 835
          }
835 836
        }
836 837

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

	
19 19
#ifndef LEMON_COST_SCALING_H
20 20
#define LEMON_COST_SCALING_H
21 21

	
22 22
/// \ingroup min_cost_flow_algs
23 23
/// \file
24 24
/// \brief Cost scaling algorithm for finding a minimum cost flow.
25 25

	
26 26
#include <vector>
27 27
#include <deque>
28 28
#include <limits>
29 29

	
30 30
#include <lemon/core.h>
31 31
#include <lemon/maps.h>
32 32
#include <lemon/math.h>
33 33
#include <lemon/static_graph.h>
34 34
#include <lemon/circulation.h>
35 35
#include <lemon/bellman_ford.h>
36 36

	
37 37
namespace lemon {
38 38

	
39 39
  /// \brief Default traits class of CostScaling algorithm.
40 40
  ///
41 41
  /// Default traits class of CostScaling algorithm.
42 42
  /// \tparam GR Digraph type.
43 43
  /// \tparam V The number type used for flow amounts, capacity bounds
44 44
  /// and supply values. By default it is \c int.
45 45
  /// \tparam C The number type used for costs and potentials.
46 46
  /// By default it is the same as \c V.
47 47
#ifdef DOXYGEN
48 48
  template <typename GR, typename V = int, typename C = V>
49 49
#else
50 50
  template < typename GR, typename V = int, typename C = V,
51 51
             bool integer = std::numeric_limits<C>::is_integer >
52 52
#endif
53 53
  struct CostScalingDefaultTraits
54 54
  {
55 55
    /// The type of the digraph
56 56
    typedef GR Digraph;
57 57
    /// The type of the flow amounts, capacity bounds and supply values
58 58
    typedef V Value;
59 59
    /// The type of the arc costs
60 60
    typedef C Cost;
61 61

	
62 62
    /// \brief The large cost type used for internal computations
63 63
    ///
64 64
    /// The large cost type used for internal computations.
65 65
    /// It is \c long \c long if the \c Cost type is integer,
66 66
    /// otherwise it is \c double.
67 67
    /// \c Cost must be convertible to \c LargeCost.
68 68
    typedef double LargeCost;
69 69
  };
70 70

	
71 71
  // Default traits class for integer cost types
72 72
  template <typename GR, typename V, typename C>
73 73
  struct CostScalingDefaultTraits<GR, V, C, true>
74 74
  {
75 75
    typedef GR Digraph;
76 76
    typedef V Value;
77 77
    typedef C Cost;
78 78
#ifdef LEMON_HAVE_LONG_LONG
79 79
    typedef long long LargeCost;
80 80
#else
81 81
    typedef long LargeCost;
82 82
#endif
83 83
  };
84 84

	
85 85

	
86 86
  /// \addtogroup min_cost_flow_algs
87 87
  /// @{
88 88

	
89 89
  /// \brief Implementation of the Cost Scaling algorithm for
90 90
  /// finding a \ref min_cost_flow "minimum cost flow".
91 91
  ///
92 92
  /// \ref CostScaling implements a cost scaling algorithm that performs
93
  /// push/augment and relabel operations for finding a minimum cost
94
  /// flow. It is an efficient primal-dual solution method, which
93
  /// push/augment and relabel operations for finding a \ref min_cost_flow
94
  /// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation,
95
  /// \ref goldberg97efficient, \ref bunnagel98efficient. 
96
  /// It is a highly efficient primal-dual solution method, which
95 97
  /// can be viewed as the generalization of the \ref Preflow
96 98
  /// "preflow push-relabel" algorithm for the maximum flow problem.
97 99
  ///
98 100
  /// Most of the parameters of the problem (except for the digraph)
99 101
  /// can be given using separate functions, and the algorithm can be
100 102
  /// executed using the \ref run() function. If some parameters are not
101 103
  /// specified, then default values will be used.
102 104
  ///
103 105
  /// \tparam GR The digraph type the algorithm runs on.
104 106
  /// \tparam V The number type used for flow amounts, capacity bounds
105 107
  /// and supply values in the algorithm. By default it is \c int.
106 108
  /// \tparam C The number type used for costs and potentials in the
107 109
  /// algorithm. By default it is the same as \c V.
108 110
  ///
109 111
  /// \warning Both number types must be signed and all input data must
110 112
  /// be integer.
111 113
  /// \warning This algorithm does not support negative costs for such
112 114
  /// arcs that have infinite upper bound.
113 115
  ///
114 116
  /// \note %CostScaling provides three different internal methods,
115 117
  /// from which the most efficient one is used by default.
116 118
  /// For more information, see \ref Method.
117 119
#ifdef DOXYGEN
118 120
  template <typename GR, typename V, typename C, typename TR>
119 121
#else
120 122
  template < typename GR, typename V = int, typename C = V,
121 123
             typename TR = CostScalingDefaultTraits<GR, V, C> >
122 124
#endif
123 125
  class CostScaling
124 126
  {
125 127
  public:
126 128

	
127 129
    /// The type of the digraph
128 130
    typedef typename TR::Digraph Digraph;
129 131
    /// The type of the flow amounts, capacity bounds and supply values
130 132
    typedef typename TR::Value Value;
131 133
    /// The type of the arc costs
132 134
    typedef typename TR::Cost Cost;
133 135

	
134 136
    /// \brief The large cost type
135 137
    ///
136 138
    /// The large cost type used for internal computations.
137 139
    /// Using the \ref CostScalingDefaultTraits "default traits class",
138 140
    /// it is \c long \c long if the \c Cost type is integer,
139 141
    /// otherwise it is \c double.
140 142
    typedef typename TR::LargeCost LargeCost;
141 143

	
142 144
    /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
143 145
    typedef TR Traits;
144 146

	
145 147
  public:
146 148

	
147 149
    /// \brief Problem type constants for the \c run() function.
148 150
    ///
149 151
    /// Enum type containing the problem type constants that can be
150 152
    /// returned by the \ref run() function of the algorithm.
151 153
    enum ProblemType {
152 154
      /// The problem has no feasible solution (flow).
153 155
      INFEASIBLE,
154 156
      /// The problem has optimal solution (i.e. it is feasible and
155 157
      /// bounded), and the algorithm has found optimal flow and node
156 158
      /// potentials (primal and dual solutions).
157 159
      OPTIMAL,
158 160
      /// The digraph contains an arc of negative cost and infinite
159 161
      /// upper bound. It means that the objective function is unbounded
160 162
      /// on that arc, however, note that it could actually be bounded
161 163
      /// over the feasible flows, but this algroithm cannot handle
162 164
      /// these cases.
163 165
      UNBOUNDED
164 166
    };
165 167

	
166 168
    /// \brief Constants for selecting the internal method.
167 169
    ///
168 170
    /// Enum type containing constants for selecting the internal method
169 171
    /// for the \ref run() function.
170 172
    ///
171 173
    /// \ref CostScaling provides three internal methods that differ mainly
172 174
    /// in their base operations, which are used in conjunction with the
173 175
    /// relabel operation.
174 176
    /// By default, the so called \ref PARTIAL_AUGMENT
175 177
    /// "Partial Augment-Relabel" method is used, which proved to be
176 178
    /// the most efficient and the most robust on various test inputs.
177 179
    /// However, the other methods can be selected using the \ref run()
178 180
    /// function with the proper parameter.
179 181
    enum Method {
180 182
      /// Local push operations are used, i.e. flow is moved only on one
181 183
      /// admissible arc at once.
182 184
      PUSH,
183 185
      /// Augment operations are used, i.e. flow is moved on admissible
184 186
      /// paths from a node with excess to a node with deficit.
185 187
      AUGMENT,
186 188
      /// Partial augment operations are used, i.e. flow is moved on 
187 189
      /// admissible paths started from a node with excess, but the
188 190
      /// lengths of these paths are limited. This method can be viewed
189 191
      /// as a combined version of the previous two operations.
190 192
      PARTIAL_AUGMENT
191 193
    };
192 194

	
193 195
  private:
194 196

	
195 197
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
196 198

	
197 199
    typedef std::vector<int> IntVector;
198 200
    typedef std::vector<char> BoolVector;
199 201
    typedef std::vector<Value> ValueVector;
200 202
    typedef std::vector<Cost> CostVector;
201 203
    typedef std::vector<LargeCost> LargeCostVector;
202 204

	
203 205
  private:
204 206
  
205 207
    template <typename KT, typename VT>
206 208
    class VectorMap {
207 209
    public:
208 210
      typedef KT Key;
209 211
      typedef VT Value;
210 212
      
211 213
      VectorMap(std::vector<Value>& v) : _v(v) {}
212 214
      
213 215
      const Value& operator[](const Key& key) const {
214 216
        return _v[StaticDigraph::id(key)];
215 217
      }
216 218

	
217 219
      Value& operator[](const Key& key) {
218 220
        return _v[StaticDigraph::id(key)];
219 221
      }
220 222
      
221 223
      void set(const Key& key, const Value& val) {
222 224
        _v[StaticDigraph::id(key)] = val;
223 225
      }
224 226

	
225 227
    private:
226 228
      std::vector<Value>& _v;
227 229
    };
228 230

	
229 231
    typedef VectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
230 232
    typedef VectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
231 233

	
232 234
  private:
233 235

	
234 236
    // Data related to the underlying digraph
235 237
    const GR &_graph;
236 238
    int _node_num;
237 239
    int _arc_num;
238 240
    int _res_node_num;
239 241
    int _res_arc_num;
240 242
    int _root;
241 243

	
242 244
    // Parameters of the problem
243 245
    bool _have_lower;
244 246
    Value _sum_supply;
245 247

	
246 248
    // Data structures for storing the digraph
247 249
    IntNodeMap _node_id;
248 250
    IntArcMap _arc_idf;
249 251
    IntArcMap _arc_idb;
250 252
    IntVector _first_out;
251 253
    BoolVector _forward;
252 254
    IntVector _source;
253 255
    IntVector _target;
254 256
    IntVector _reverse;
255 257

	
256 258
    // Node and arc data
257 259
    ValueVector _lower;
258 260
    ValueVector _upper;
259 261
    CostVector _scost;
260 262
    ValueVector _supply;
261 263

	
262 264
    ValueVector _res_cap;
263 265
    LargeCostVector _cost;
264 266
    LargeCostVector _pi;
265 267
    ValueVector _excess;
266 268
    IntVector _next_out;
267 269
    std::deque<int> _active_nodes;
268 270

	
269 271
    // Data for scaling
270 272
    LargeCost _epsilon;
271 273
    int _alpha;
272 274

	
273 275
    // Data for a StaticDigraph structure
274 276
    typedef std::pair<int, int> IntPair;
275 277
    StaticDigraph _sgr;
276 278
    std::vector<IntPair> _arc_vec;
277 279
    std::vector<LargeCost> _cost_vec;
278 280
    LargeCostArcMap _cost_map;
279 281
    LargeCostNodeMap _pi_map;
280 282
  
281 283
  public:
282 284
  
283 285
    /// \brief Constant for infinite upper bounds (capacities).
284 286
    ///
285 287
    /// Constant for infinite upper bounds (capacities).
286 288
    /// It is \c std::numeric_limits<Value>::infinity() if available,
287 289
    /// \c std::numeric_limits<Value>::max() otherwise.
288 290
    const Value INF;
289 291

	
290 292
  public:
291 293

	
292 294
    /// \name Named Template Parameters
293 295
    /// @{
294 296

	
295 297
    template <typename T>
296 298
    struct SetLargeCostTraits : public Traits {
297 299
      typedef T LargeCost;
298 300
    };
299 301

	
300 302
    /// \brief \ref named-templ-param "Named parameter" for setting
301 303
    /// \c LargeCost type.
302 304
    ///
303 305
    /// \ref named-templ-param "Named parameter" for setting \c LargeCost
304 306
    /// type, which is used for internal computations in the algorithm.
305 307
    /// \c Cost must be convertible to \c LargeCost.
306 308
    template <typename T>
307 309
    struct SetLargeCost
308 310
      : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
309 311
      typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
310 312
    };
311 313

	
312 314
    /// @}
313 315

	
314 316
  public:
315 317

	
316 318
    /// \brief Constructor.
317 319
    ///
318 320
    /// The constructor of the class.
319 321
    ///
320 322
    /// \param graph The digraph the algorithm runs on.
321 323
    CostScaling(const GR& graph) :
322 324
      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
323 325
      _cost_map(_cost_vec), _pi_map(_pi),
324 326
      INF(std::numeric_limits<Value>::has_infinity ?
325 327
          std::numeric_limits<Value>::infinity() :
326 328
          std::numeric_limits<Value>::max())
327 329
    {
328 330
      // Check the number types
329 331
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
330 332
        "The flow type of CostScaling must be signed");
331 333
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
332 334
        "The cost type of CostScaling must be signed");
333 335

	
334 336
      // Resize vectors
335 337
      _node_num = countNodes(_graph);
336 338
      _arc_num = countArcs(_graph);
337 339
      _res_node_num = _node_num + 1;
338 340
      _res_arc_num = 2 * (_arc_num + _node_num);
339 341
      _root = _node_num;
340 342

	
341 343
      _first_out.resize(_res_node_num + 1);
342 344
      _forward.resize(_res_arc_num);
343 345
      _source.resize(_res_arc_num);
344 346
      _target.resize(_res_arc_num);
345 347
      _reverse.resize(_res_arc_num);
346 348

	
347 349
      _lower.resize(_res_arc_num);
348 350
      _upper.resize(_res_arc_num);
349 351
      _scost.resize(_res_arc_num);
350 352
      _supply.resize(_res_node_num);
351 353
      
352 354
      _res_cap.resize(_res_arc_num);
353 355
      _cost.resize(_res_arc_num);
354 356
      _pi.resize(_res_node_num);
355 357
      _excess.resize(_res_node_num);
356 358
      _next_out.resize(_res_node_num);
357 359

	
358 360
      _arc_vec.reserve(_res_arc_num);
359 361
      _cost_vec.reserve(_res_arc_num);
360 362

	
361 363
      // Copy the graph
362 364
      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
363 365
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
364 366
        _node_id[n] = i;
365 367
      }
366 368
      i = 0;
367 369
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
368 370
        _first_out[i] = j;
369 371
        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
370 372
          _arc_idf[a] = j;
371 373
          _forward[j] = true;
372 374
          _source[j] = i;
373 375
          _target[j] = _node_id[_graph.runningNode(a)];
374 376
        }
375 377
        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
376 378
          _arc_idb[a] = j;
377 379
          _forward[j] = false;
378 380
          _source[j] = i;
379 381
          _target[j] = _node_id[_graph.runningNode(a)];
380 382
        }
381 383
        _forward[j] = false;
382 384
        _source[j] = i;
383 385
        _target[j] = _root;
384 386
        _reverse[j] = k;
385 387
        _forward[k] = true;
386 388
        _source[k] = _root;
387 389
        _target[k] = i;
388 390
        _reverse[k] = j;
389 391
        ++j; ++k;
390 392
      }
391 393
      _first_out[i] = j;
392 394
      _first_out[_res_node_num] = k;
393 395
      for (ArcIt a(_graph); a != INVALID; ++a) {
394 396
        int fi = _arc_idf[a];
395 397
        int bi = _arc_idb[a];
396 398
        _reverse[fi] = bi;
397 399
        _reverse[bi] = fi;
398 400
      }
399 401
      
400 402
      // Reset parameters
401 403
      reset();
402 404
    }
403 405

	
404 406
    /// \name Parameters
405 407
    /// The parameters of the algorithm can be specified using these
406 408
    /// functions.
407 409

	
408 410
    /// @{
409 411

	
410 412
    /// \brief Set the lower bounds on the arcs.
411 413
    ///
412 414
    /// This function sets the lower bounds on the arcs.
413 415
    /// If it is not used before calling \ref run(), the lower bounds
414 416
    /// will be set to zero on all arcs.
415 417
    ///
416 418
    /// \param map An arc map storing the lower bounds.
417 419
    /// Its \c Value type must be convertible to the \c Value type
418 420
    /// of the algorithm.
419 421
    ///
420 422
    /// \return <tt>(*this)</tt>
421 423
    template <typename LowerMap>
422 424
    CostScaling& lowerMap(const LowerMap& map) {
423 425
      _have_lower = true;
424 426
      for (ArcIt a(_graph); a != INVALID; ++a) {
425 427
        _lower[_arc_idf[a]] = map[a];
426 428
        _lower[_arc_idb[a]] = map[a];
427 429
      }
428 430
      return *this;
429 431
    }
430 432

	
431 433
    /// \brief Set the upper bounds (capacities) on the arcs.
432 434
    ///
433 435
    /// This function sets the upper bounds (capacities) on the arcs.
434 436
    /// If it is not used before calling \ref run(), the upper bounds
435 437
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
436 438
    /// unbounded from above).
437 439
    ///
438 440
    /// \param map An arc map storing the upper bounds.
439 441
    /// Its \c Value type must be convertible to the \c Value type
440 442
    /// of the algorithm.
441 443
    ///
442 444
    /// \return <tt>(*this)</tt>
443 445
    template<typename UpperMap>
444 446
    CostScaling& upperMap(const UpperMap& map) {
445 447
      for (ArcIt a(_graph); a != INVALID; ++a) {
446 448
        _upper[_arc_idf[a]] = map[a];
447 449
      }
448 450
      return *this;
449 451
    }
450 452

	
451 453
    /// \brief Set the costs of the arcs.
452 454
    ///
453 455
    /// This function sets the costs of the arcs.
454 456
    /// If it is not used before calling \ref run(), the costs
455 457
    /// will be set to \c 1 on all arcs.
456 458
    ///
457 459
    /// \param map An arc map storing the costs.
458 460
    /// Its \c Value type must be convertible to the \c Cost type
459 461
    /// of the algorithm.
460 462
    ///
461 463
    /// \return <tt>(*this)</tt>
462 464
    template<typename CostMap>
463 465
    CostScaling& costMap(const CostMap& map) {
464 466
      for (ArcIt a(_graph); a != INVALID; ++a) {
465 467
        _scost[_arc_idf[a]] =  map[a];
466 468
        _scost[_arc_idb[a]] = -map[a];
467 469
      }
468 470
      return *this;
469 471
    }
470 472

	
471 473
    /// \brief Set the supply values of the nodes.
472 474
    ///
473 475
    /// This function sets the supply values of the nodes.
474 476
    /// If neither this function nor \ref stSupply() is used before
475 477
    /// calling \ref run(), the supply of each node will be set to zero.
476 478
    ///
477 479
    /// \param map A node map storing the supply values.
478 480
    /// Its \c Value type must be convertible to the \c Value type
479 481
    /// of the algorithm.
480 482
    ///
481 483
    /// \return <tt>(*this)</tt>
482 484
    template<typename SupplyMap>
483 485
    CostScaling& supplyMap(const SupplyMap& map) {
484 486
      for (NodeIt n(_graph); n != INVALID; ++n) {
485 487
        _supply[_node_id[n]] = map[n];
486 488
      }
487 489
      return *this;
488 490
    }
489 491

	
490 492
    /// \brief Set single source and target nodes and a supply value.
491 493
    ///
492 494
    /// This function sets a single source node and a single target node
493 495
    /// and the required flow value.
494 496
    /// If neither this function nor \ref supplyMap() is used before
495 497
    /// calling \ref run(), the supply of each node will be set to zero.
496 498
    ///
497 499
    /// Using this function has the same effect as using \ref supplyMap()
498 500
    /// with such a map in which \c k is assigned to \c s, \c -k is
499 501
    /// assigned to \c t and all other nodes have zero supply value.
500 502
    ///
501 503
    /// \param s The source node.
502 504
    /// \param t The target node.
503 505
    /// \param k The required amount of flow from node \c s to node \c t
504 506
    /// (i.e. the supply of \c s and the demand of \c t).
505 507
    ///
506 508
    /// \return <tt>(*this)</tt>
507 509
    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
508 510
      for (int i = 0; i != _res_node_num; ++i) {
509 511
        _supply[i] = 0;
510 512
      }
511 513
      _supply[_node_id[s]] =  k;
512 514
      _supply[_node_id[t]] = -k;
513 515
      return *this;
514 516
    }
515 517
    
516 518
    /// @}
517 519

	
518 520
    /// \name Execution control
519 521
    /// The algorithm can be executed using \ref run().
520 522

	
521 523
    /// @{
522 524

	
523 525
    /// \brief Run the algorithm.
524 526
    ///
525 527
    /// This function runs the algorithm.
526 528
    /// The paramters can be specified using functions \ref lowerMap(),
527 529
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
528 530
    /// For example,
529 531
    /// \code
530 532
    ///   CostScaling<ListDigraph> cs(graph);
531 533
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
532 534
    ///     .supplyMap(sup).run();
533 535
    /// \endcode
534 536
    ///
535 537
    /// This function can be called more than once. All the parameters
536 538
    /// that have been given are kept for the next call, unless
537 539
    /// \ref reset() is called, thus only the modified parameters
538 540
    /// have to be set again. See \ref reset() for examples.
539 541
    /// However, the underlying digraph must not be modified after this
540 542
    /// class have been constructed, since it copies and extends the graph.
541 543
    ///
542 544
    /// \param method The internal method that will be used in the
543 545
    /// algorithm. For more information, see \ref Method.
544 546
    /// \param factor The cost scaling factor. It must be larger than one.
545 547
    ///
546 548
    /// \return \c INFEASIBLE if no feasible flow exists,
547 549
    /// \n \c OPTIMAL if the problem has optimal solution
548 550
    /// (i.e. it is feasible and bounded), and the algorithm has found
549 551
    /// optimal flow and node potentials (primal and dual solutions),
550 552
    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
551 553
    /// and infinite upper bound. It means that the objective function
552 554
    /// is unbounded on that arc, however, note that it could actually be
553 555
    /// bounded over the feasible flows, but this algroithm cannot handle
554 556
    /// these cases.
555 557
    ///
556 558
    /// \see ProblemType, Method
557 559
    ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
558 560
      _alpha = factor;
559 561
      ProblemType pt = init();
560 562
      if (pt != OPTIMAL) return pt;
561 563
      start(method);
562 564
      return OPTIMAL;
563 565
    }
564 566

	
565 567
    /// \brief Reset all the parameters that have been given before.
566 568
    ///
567 569
    /// This function resets all the paramaters that have been given
568 570
    /// before using functions \ref lowerMap(), \ref upperMap(),
569 571
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
570 572
    ///
571 573
    /// It is useful for multiple run() calls. If this function is not
572 574
    /// used, all the parameters given before are kept for the next
573 575
    /// \ref run() call.
574 576
    /// However, the underlying digraph must not be modified after this
575 577
    /// class have been constructed, since it copies and extends the graph.
576 578
    ///
577 579
    /// For example,
578 580
    /// \code
579 581
    ///   CostScaling<ListDigraph> cs(graph);
580 582
    ///
581 583
    ///   // First run
582 584
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
583 585
    ///     .supplyMap(sup).run();
584 586
    ///
585 587
    ///   // Run again with modified cost map (reset() is not called,
586 588
    ///   // so only the cost map have to be set again)
587 589
    ///   cost[e] += 100;
588 590
    ///   cs.costMap(cost).run();
589 591
    ///
590 592
    ///   // Run again from scratch using reset()
591 593
    ///   // (the lower bounds will be set to zero on all arcs)
592 594
    ///   cs.reset();
593 595
    ///   cs.upperMap(capacity).costMap(cost)
594 596
    ///     .supplyMap(sup).run();
595 597
    /// \endcode
596 598
    ///
597 599
    /// \return <tt>(*this)</tt>
598 600
    CostScaling& reset() {
599 601
      for (int i = 0; i != _res_node_num; ++i) {
600 602
        _supply[i] = 0;
601 603
      }
602 604
      int limit = _first_out[_root];
603 605
      for (int j = 0; j != limit; ++j) {
604 606
        _lower[j] = 0;
605 607
        _upper[j] = INF;
606 608
        _scost[j] = _forward[j] ? 1 : -1;
607 609
      }
608 610
      for (int j = limit; j != _res_arc_num; ++j) {
609 611
        _lower[j] = 0;
610 612
        _upper[j] = INF;
611 613
        _scost[j] = 0;
612 614
        _scost[_reverse[j]] = 0;
613 615
      }      
614 616
      _have_lower = false;
615 617
      return *this;
616 618
    }
617 619

	
618 620
    /// @}
619 621

	
620 622
    /// \name Query Functions
621 623
    /// The results of the algorithm can be obtained using these
622 624
    /// functions.\n
623 625
    /// The \ref run() function must be called before using them.
624 626

	
625 627
    /// @{
626 628

	
627 629
    /// \brief Return the total cost of the found flow.
628 630
    ///
629 631
    /// This function returns the total cost of the found flow.
630 632
    /// Its complexity is O(e).
631 633
    ///
632 634
    /// \note The return type of the function can be specified as a
633 635
    /// template parameter. For example,
634 636
    /// \code
635 637
    ///   cs.totalCost<double>();
636 638
    /// \endcode
637 639
    /// It is useful if the total cost cannot be stored in the \c Cost
638 640
    /// type of the algorithm, which is the default return type of the
639 641
    /// function.
640 642
    ///
641 643
    /// \pre \ref run() must be called before using this function.
642 644
    template <typename Number>
643 645
    Number totalCost() const {
644 646
      Number c = 0;
645 647
      for (ArcIt a(_graph); a != INVALID; ++a) {
646 648
        int i = _arc_idb[a];
647 649
        c += static_cast<Number>(_res_cap[i]) *
648 650
             (-static_cast<Number>(_scost[i]));
649 651
      }
650 652
      return c;
651 653
    }
652 654

	
653 655
#ifndef DOXYGEN
654 656
    Cost totalCost() const {
655 657
      return totalCost<Cost>();
656 658
    }
657 659
#endif
658 660

	
659 661
    /// \brief Return the flow on the given arc.
660 662
    ///
661 663
    /// This function returns the flow on the given arc.
662 664
    ///
663 665
    /// \pre \ref run() must be called before using this function.
664 666
    Value flow(const Arc& a) const {
665 667
      return _res_cap[_arc_idb[a]];
666 668
    }
667 669

	
668 670
    /// \brief Return the flow map (the primal solution).
669 671
    ///
670 672
    /// This function copies the flow value on each arc into the given
671 673
    /// map. The \c Value type of the algorithm must be convertible to
672 674
    /// the \c Value type of the map.
673 675
    ///
674 676
    /// \pre \ref run() must be called before using this function.
675 677
    template <typename FlowMap>
676 678
    void flowMap(FlowMap &map) const {
677 679
      for (ArcIt a(_graph); a != INVALID; ++a) {
678 680
        map.set(a, _res_cap[_arc_idb[a]]);
679 681
      }
680 682
    }
681 683

	
682 684
    /// \brief Return the potential (dual value) of the given node.
683 685
    ///
684 686
    /// This function returns the potential (dual value) of the
685 687
    /// given node.
686 688
    ///
687 689
    /// \pre \ref run() must be called before using this function.
688 690
    Cost potential(const Node& n) const {
689 691
      return static_cast<Cost>(_pi[_node_id[n]]);
690 692
    }
691 693

	
692 694
    /// \brief Return the potential map (the dual solution).
693 695
    ///
694 696
    /// This function copies the potential (dual value) of each node
695 697
    /// into the given map.
696 698
    /// The \c Cost type of the algorithm must be convertible to the
697 699
    /// \c Value type of the map.
698 700
    ///
699 701
    /// \pre \ref run() must be called before using this function.
700 702
    template <typename PotentialMap>
701 703
    void potentialMap(PotentialMap &map) const {
702 704
      for (NodeIt n(_graph); n != INVALID; ++n) {
703 705
        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
704 706
      }
705 707
    }
706 708

	
707 709
    /// @}
708 710

	
709 711
  private:
710 712

	
711 713
    // Initialize the algorithm
712 714
    ProblemType init() {
713 715
      if (_res_node_num == 0) return INFEASIBLE;
714 716

	
715 717
      // Check the sum of supply values
716 718
      _sum_supply = 0;
717 719
      for (int i = 0; i != _root; ++i) {
718 720
        _sum_supply += _supply[i];
719 721
      }
720 722
      if (_sum_supply > 0) return INFEASIBLE;
721 723
      
722 724

	
723 725
      // Initialize vectors
724 726
      for (int i = 0; i != _res_node_num; ++i) {
725 727
        _pi[i] = 0;
726 728
        _excess[i] = _supply[i];
727 729
      }
728 730
      
729 731
      // Remove infinite upper bounds and check negative arcs
730 732
      const Value MAX = std::numeric_limits<Value>::max();
731 733
      int last_out;
732 734
      if (_have_lower) {
733 735
        for (int i = 0; i != _root; ++i) {
734 736
          last_out = _first_out[i+1];
735 737
          for (int j = _first_out[i]; j != last_out; ++j) {
736 738
            if (_forward[j]) {
737 739
              Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
738 740
              if (c >= MAX) return UNBOUNDED;
739 741
              _excess[i] -= c;
740 742
              _excess[_target[j]] += c;
741 743
            }
742 744
          }
743 745
        }
744 746
      } else {
745 747
        for (int i = 0; i != _root; ++i) {
746 748
          last_out = _first_out[i+1];
747 749
          for (int j = _first_out[i]; j != last_out; ++j) {
748 750
            if (_forward[j] && _scost[j] < 0) {
749 751
              Value c = _upper[j];
750 752
              if (c >= MAX) return UNBOUNDED;
751 753
              _excess[i] -= c;
752 754
              _excess[_target[j]] += c;
753 755
            }
754 756
          }
755 757
        }
756 758
      }
757 759
      Value ex, max_cap = 0;
758 760
      for (int i = 0; i != _res_node_num; ++i) {
759 761
        ex = _excess[i];
760 762
        _excess[i] = 0;
761 763
        if (ex < 0) max_cap -= ex;
762 764
      }
763 765
      for (int j = 0; j != _res_arc_num; ++j) {
764 766
        if (_upper[j] >= MAX) _upper[j] = max_cap;
765 767
      }
766 768

	
767 769
      // Initialize the large cost vector and the epsilon parameter
768 770
      _epsilon = 0;
769 771
      LargeCost lc;
770 772
      for (int i = 0; i != _root; ++i) {
771 773
        last_out = _first_out[i+1];
772 774
        for (int j = _first_out[i]; j != last_out; ++j) {
773 775
          lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
774 776
          _cost[j] = lc;
775 777
          if (lc > _epsilon) _epsilon = lc;
776 778
        }
777 779
      }
778 780
      _epsilon /= _alpha;
779 781

	
780 782
      // Initialize maps for Circulation and remove non-zero lower bounds
781 783
      ConstMap<Arc, Value> low(0);
782 784
      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
783 785
      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
784 786
      ValueArcMap cap(_graph), flow(_graph);
785 787
      ValueNodeMap sup(_graph);
786 788
      for (NodeIt n(_graph); n != INVALID; ++n) {
787 789
        sup[n] = _supply[_node_id[n]];
788 790
      }
789 791
      if (_have_lower) {
790 792
        for (ArcIt a(_graph); a != INVALID; ++a) {
791 793
          int j = _arc_idf[a];
792 794
          Value c = _lower[j];
793 795
          cap[a] = _upper[j] - c;
794 796
          sup[_graph.source(a)] -= c;
795 797
          sup[_graph.target(a)] += c;
796 798
        }
797 799
      } else {
798 800
        for (ArcIt a(_graph); a != INVALID; ++a) {
799 801
          cap[a] = _upper[_arc_idf[a]];
800 802
        }
801 803
      }
802 804

	
803 805
      // Find a feasible flow using Circulation
804 806
      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
805 807
        circ(_graph, low, cap, sup);
806 808
      if (!circ.flowMap(flow).run()) return INFEASIBLE;
807 809

	
808 810
      // Set residual capacities and handle GEQ supply type
809 811
      if (_sum_supply < 0) {
810 812
        for (ArcIt a(_graph); a != INVALID; ++a) {
811 813
          Value fa = flow[a];
812 814
          _res_cap[_arc_idf[a]] = cap[a] - fa;
813 815
          _res_cap[_arc_idb[a]] = fa;
814 816
          sup[_graph.source(a)] -= fa;
815 817
          sup[_graph.target(a)] += fa;
816 818
        }
817 819
        for (NodeIt n(_graph); n != INVALID; ++n) {
818 820
          _excess[_node_id[n]] = sup[n];
819 821
        }
820 822
        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
821 823
          int u = _target[a];
822 824
          int ra = _reverse[a];
823 825
          _res_cap[a] = -_sum_supply + 1;
824 826
          _res_cap[ra] = -_excess[u];
825 827
          _cost[a] = 0;
826 828
          _cost[ra] = 0;
827 829
          _excess[u] = 0;
828 830
        }
829 831
      } else {
830 832
        for (ArcIt a(_graph); a != INVALID; ++a) {
831 833
          Value fa = flow[a];
832 834
          _res_cap[_arc_idf[a]] = cap[a] - fa;
833 835
          _res_cap[_arc_idb[a]] = fa;
834 836
        }
835 837
        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
836 838
          int ra = _reverse[a];
837 839
          _res_cap[a] = 1;
838 840
          _res_cap[ra] = 0;
839 841
          _cost[a] = 0;
840 842
          _cost[ra] = 0;
841 843
        }
842 844
      }
843 845
      
844 846
      return OPTIMAL;
845 847
    }
846 848

	
847 849
    // Execute the algorithm and transform the results
848 850
    void start(Method method) {
849 851
      // Maximum path length for partial augment
850 852
      const int MAX_PATH_LENGTH = 4;
851 853
      
852 854
      // Execute the algorithm
853 855
      switch (method) {
854 856
        case PUSH:
855 857
          startPush();
856 858
          break;
857 859
        case AUGMENT:
858 860
          startAugment();
859 861
          break;
860 862
        case PARTIAL_AUGMENT:
861 863
          startAugment(MAX_PATH_LENGTH);
862 864
          break;
0 comments (0 inline)