gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Various doc improvements (#406)
0 8 0
default
8 files changed with 32 insertions and 28 deletions:
↑ Collapse diff ↑
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
/*!
20 20

	
21 21
\page coding_style LEMON Coding Style
22 22

	
23 23
\section naming_conv Naming Conventions
24 24

	
25 25
In order to make development easier we have made some conventions
26 26
according to coding style. These include names of types, classes,
27 27
functions, variables, constants and exceptions. If these conventions
28 28
are met in one's code then it is easier to read and maintain
29 29
it. Please comply with these conventions if you want to contribute
30 30
developing LEMON library.
31 31

	
32 32
\note When the coding style requires the capitalization of an abbreviation,
33 33
only the first letter should be upper case.
34 34

	
35 35
\code
36 36
XmlReader
37 37
\endcode
38 38

	
39 39

	
40 40
\warning In some cases we diverge from these rules.
41 41
This is primary done because STL uses different naming convention and
42 42
in certain cases
43 43
it is beneficial to provide STL compatible interface.
44 44

	
45 45
\subsection cs-files File Names
46 46

	
47 47
The header file names should look like the following.
48 48

	
49 49
\code
50 50
header_file.h
51 51
\endcode
52 52

	
53 53
Note that all standard LEMON headers are located in the \c lemon subdirectory,
54 54
so you should include them from C++ source like this:
55 55

	
56 56
\code
57 57
#include <lemon/header_file.h>
58 58
\endcode
59 59

	
60 60
The source code files use the same style and they have '.cc' extension.
61 61

	
62 62
\code
63 63
source_code.cc
64 64
\endcode
65 65

	
66 66
\subsection cs-class Classes and other types
67 67

	
68 68
The name of a class or any type should look like the following.
69 69

	
70 70
\code
71 71
AllWordsCapitalizedWithoutUnderscores
72 72
\endcode
73 73

	
74 74
\subsection cs-func Methods and other functions
75 75

	
76 76
The name of a function should look like the following.
77 77

	
78 78
\code
79 79
firstWordLowerCaseRestCapitalizedWithoutUnderscores
80 80
\endcode
81 81

	
82 82
\subsection cs-funcs Constants, Macros
83 83

	
84 84
The names of constants and macros should look like the following.
85 85

	
86 86
\code
87 87
ALL_UPPER_CASE_WITH_UNDERSCORES
88 88
\endcode
89 89

	
90 90
\subsection cs-loc-var Class and instance member variables, auto variables
91 91

	
92 92
The names of class and instance member variables and auto variables
93 93
(=variables used locally in methods) should look like the following.
94 94

	
95 95
\code
96 96
all_lower_case_with_underscores
97 97
\endcode
98 98

	
99 99
\subsection pri-loc-var Private member variables
100 100

	
101
Private member variables should start with underscore
101
Private member variables should start with underscore.
102 102

	
103 103
\code
104
_start_with_underscores
104
_start_with_underscore
105 105
\endcode
106 106

	
107 107
\subsection cs-excep Exceptions
108 108

	
109 109
When writing exceptions please comply the following naming conventions.
110 110

	
111 111
\code
112 112
ClassNameEndsWithException
113 113
\endcode
114 114

	
115 115
or
116 116

	
117 117
\code
118 118
ClassNameEndsWithError
119 119
\endcode
120 120

	
121 121
\section header-template Template Header File
122 122

	
123 123
Each LEMON header file should look like this:
124 124

	
125 125
\include template.h
126 126

	
127 127
*/
Ignore white space 6 line context
... ...
@@ -217,518 +217,518 @@
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 auxdat Auxiliary Data Structures
267 267
@ingroup datas
268 268
\brief Auxiliary data structures implemented in LEMON.
269 269

	
270 270
This group contains some data structures implemented in LEMON in
271 271
order to make it easier to implement combinatorial algorithms.
272 272
*/
273 273

	
274 274
/**
275 275
@defgroup geomdat Geometric Data Structures
276 276
@ingroup auxdat
277 277
\brief Geometric data structures implemented in LEMON.
278 278

	
279 279
This group contains geometric data structures implemented in LEMON.
280 280

	
281 281
 - \ref lemon::dim2::Point "dim2::Point" implements a two dimensional
282 282
   vector with the usual operations.
283 283
 - \ref lemon::dim2::Box "dim2::Box" can be used to determine the
284 284
   rectangular bounding box of a set of \ref lemon::dim2::Point
285 285
   "dim2::Point"'s.
286 286
*/
287 287

	
288 288
/**
289 289
@defgroup matrices Matrices
290 290
@ingroup auxdat
291 291
\brief Two dimensional data storages implemented in LEMON.
292 292

	
293 293
This group contains two dimensional data storages implemented in LEMON.
294 294
*/
295 295

	
296 296
/**
297 297
@defgroup algs Algorithms
298 298
\brief This group contains the several algorithms
299 299
implemented in LEMON.
300 300

	
301 301
This group contains the several algorithms
302 302
implemented in LEMON.
303 303
*/
304 304

	
305 305
/**
306 306
@defgroup search Graph Search
307 307
@ingroup algs
308 308
\brief Common graph search algorithms.
309 309

	
310 310
This group contains the common graph search algorithms, namely
311 311
\e breadth-first \e search (BFS) and \e depth-first \e search (DFS)
312 312
\ref clrs01algorithms.
313 313
*/
314 314

	
315 315
/**
316 316
@defgroup shortest_path Shortest Path Algorithms
317 317
@ingroup algs
318 318
\brief Algorithms for finding shortest paths.
319 319

	
320 320
This group contains the algorithms for finding shortest paths in digraphs
321 321
\ref clrs01algorithms.
322 322

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

	
337 337
/**
338 338
@defgroup spantree Minimum Spanning Tree Algorithms
339 339
@ingroup algs
340 340
\brief Algorithms for finding minimum cost spanning trees and arborescences.
341 341

	
342 342
This group contains the algorithms for finding minimum cost spanning
343 343
trees and arborescences \ref clrs01algorithms.
344 344
*/
345 345

	
346 346
/**
347 347
@defgroup max_flow Maximum Flow Algorithms
348 348
@ingroup algs
349 349
\brief Algorithms for finding maximum flows.
350 350

	
351 351
This group contains the algorithms for finding maximum flows and
352 352
feasible circulations \ref clrs01algorithms, \ref amo93networkflows.
353 353

	
354 354
The \e maximum \e flow \e problem is to find a flow of maximum value between
355 355
a single source and a single target. Formally, there is a \f$G=(V,A)\f$
356 356
digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and
357 357
\f$s, t \in V\f$ source and target nodes.
358 358
A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the
359 359
following optimization problem.
360 360

	
361 361
\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f]
362 362
\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu)
363 363
    \quad \forall u\in V\setminus\{s,t\} \f]
364 364
\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f]
365 365

	
366 366
LEMON contains several algorithms for solving maximum flow problems:
367 367
- \ref EdmondsKarp Edmonds-Karp algorithm
368 368
  \ref edmondskarp72theoretical.
369 369
- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm
370 370
  \ref goldberg88newapproach.
371 371
- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees
372 372
  \ref dinic70algorithm, \ref sleator83dynamic.
373 373
- \ref GoldbergTarjan !Preflow push-relabel algorithm with dynamic trees
374 374
  \ref goldberg88newapproach, \ref sleator83dynamic.
375 375

	
376 376
In most cases the \ref Preflow algorithm provides the
377 377
fastest method for computing a maximum flow. All implementations
378 378
also provide functions to query the minimum cut, which is the dual
379 379
problem of maximum flow.
380 380

	
381 381
\ref Circulation is a preflow push-relabel algorithm implemented directly
382 382
for finding feasible circulations, which is a somewhat different problem,
383 383
but it is strongly related to maximum flow.
384 384
For more information, see \ref Circulation.
385 385
*/
386 386

	
387 387
/**
388 388
@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms
389 389
@ingroup algs
390 390

	
391 391
\brief Algorithms for finding minimum cost flows and circulations.
392 392

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

	
398 398
LEMON contains several algorithms for this problem.
399 399
 - \ref NetworkSimplex Primal Network Simplex algorithm with various
400 400
   pivot strategies \ref dantzig63linearprog, \ref kellyoneill91netsimplex.
401 401
 - \ref CostScaling Cost Scaling algorithm based on push/augment and
402 402
   relabel operations \ref goldberg90approximation, \ref goldberg97efficient,
403 403
   \ref bunnagel98efficient.
404 404
 - \ref CapacityScaling Capacity Scaling algorithm based on the successive
405 405
   shortest path method \ref edmondskarp72theoretical.
406 406
 - \ref CycleCanceling Cycle-Canceling algorithms, two of which are
407 407
   strongly polynomial \ref klein67primal, \ref goldberg89cyclecanceling.
408 408

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

	
415 415
/**
416 416
@defgroup min_cut Minimum Cut Algorithms
417 417
@ingroup algs
418 418

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

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

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

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

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

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

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

	
445 445
/**
446 446
@defgroup min_mean_cycle Minimum Mean Cycle Algorithms
447 447
@ingroup algs
448 448
\brief Algorithms for finding minimum mean cycles.
449 449

	
450 450
This group contains the algorithms for finding minimum mean cycles
451 451
\ref clrs01algorithms, \ref amo93networkflows.
452 452

	
453 453
The \e minimum \e mean \e cycle \e problem is to find a directed cycle
454 454
of minimum mean length (cost) in a digraph.
455 455
The mean length of a cycle is the average length of its arcs, i.e. the
456 456
ratio between the total length of the cycle and the number of arcs on it.
457 457

	
458 458
This problem has an important connection to \e conservative \e length
459 459
\e functions, too. A length function on the arcs of a digraph is called
460 460
conservative if and only if there is no directed cycle of negative total
461 461
length. For an arbitrary length function, the negative of the minimum
462 462
cycle mean is the smallest \f$\epsilon\f$ value so that increasing the
463 463
arc lengths uniformly by \f$\epsilon\f$ results in a conservative length
464 464
function.
465 465

	
466 466
LEMON contains three algorithms for solving the minimum mean cycle problem:
467 467
- \ref KarpMmc Karp's original algorithm \ref amo93networkflows,
468 468
  \ref dasdan98minmeancycle.
469 469
- \ref HartmannOrlinMmc Hartmann-Orlin's algorithm, which is an improved
470 470
  version of Karp's algorithm \ref dasdan98minmeancycle.
471 471
- \ref HowardMmc Howard's policy iteration algorithm
472 472
  \ref dasdan98minmeancycle.
473 473

	
474
In practice, the \ref HowardMmc "Howard" algorithm proved to be by far the
474
In practice, the \ref HowardMmc "Howard" algorithm turned out to be by far the
475 475
most efficient one, though the best known theoretical bound on its running
476 476
time is exponential.
477 477
Both \ref KarpMmc "Karp" and \ref HartmannOrlinMmc "Hartmann-Orlin" algorithms
478 478
run in time O(ne) and use space O(n<sup>2</sup>+e), but the latter one is
479 479
typically faster due to the applied early termination scheme.
480 480
*/
481 481

	
482 482
/**
483 483
@defgroup matching Matching Algorithms
484 484
@ingroup algs
485 485
\brief Algorithms for finding matchings in graphs and bipartite graphs.
486 486

	
487 487
This group contains the algorithms for calculating
488 488
matchings in graphs and bipartite graphs. The general matching problem is
489 489
finding a subset of the edges for which each node has at most one incident
490 490
edge.
491 491

	
492 492
There are several different algorithms for calculate matchings in
493 493
graphs.  The matching problems in bipartite graphs are generally
494 494
easier than in general graphs. The goal of the matching optimization
495 495
can be finding maximum cardinality, maximum weight or minimum cost
496 496
matching. The search can be constrained to find perfect or
497 497
maximum cardinality matching.
498 498

	
499 499
The matching algorithms implemented in LEMON:
500 500
- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm
501 501
  for calculating maximum cardinality matching in bipartite graphs.
502 502
- \ref PrBipartiteMatching Push-relabel algorithm
503 503
  for calculating maximum cardinality matching in bipartite graphs.
504 504
- \ref MaxWeightedBipartiteMatching
505 505
  Successive shortest path algorithm for calculating maximum weighted
506 506
  matching and maximum weighted bipartite matching in bipartite graphs.
507 507
- \ref MinCostMaxBipartiteMatching
508 508
  Successive shortest path algorithm for calculating minimum cost maximum
509 509
  matching in bipartite graphs.
510 510
- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating
511 511
  maximum cardinality matching in general graphs.
512 512
- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating
513 513
  maximum weighted matching in general graphs.
514 514
- \ref MaxWeightedPerfectMatching
515 515
  Edmond's blossom shrinking algorithm for calculating maximum weighted
516 516
  perfect matching in general graphs.
517 517
- \ref MaxFractionalMatching Push-relabel algorithm for calculating
518 518
  maximum cardinality fractional matching in general graphs.
519 519
- \ref MaxWeightedFractionalMatching Augmenting path algorithm for calculating
520 520
  maximum weighted fractional matching in general graphs.
521 521
- \ref MaxWeightedPerfectFractionalMatching
522 522
  Augmenting path algorithm for calculating maximum weighted
523 523
  perfect fractional matching in general graphs.
524 524

	
525 525
\image html matching.png
526 526
\image latex matching.eps "Min Cost Perfect Matching" width=\textwidth
527 527
*/
528 528

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

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

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

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

	
546 546
This group contains the algorithms for planarity checking,
547 547
embedding and drawing.
548 548

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

	
553 553
/**
554 554
@defgroup approx_algs Approximation Algorithms
555 555
@ingroup algs
556 556
\brief Approximation algorithms.
557 557

	
558 558
This group contains the approximation and heuristic algorithms
559 559
implemented in LEMON.
560 560

	
561 561
<b>Maximum Clique Problem</b>
562 562
  - \ref GrossoLocatelliPullanMc An efficient heuristic algorithm of
563 563
    Grosso, Locatelli, and Pullan.
564 564
*/
565 565

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

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

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

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

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

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

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

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

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

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

	
611 611
This group contains some metaheuristic optimization tools.
612 612
*/
613 613

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

	
618 618
Tools and utilities for programming in LEMON.
619 619
*/
620 620

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

	
626 626
This group contains some simple basic graph utilities.
627 627
*/
628 628

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

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

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

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

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

	
652 652
This group contains the exceptions defined in LEMON.
653 653
*/
654 654

	
655 655
/**
656 656
@defgroup io_group Input-Output
657 657
\brief Graph Input-Output methods
658 658

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

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

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

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

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

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

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

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

	
696 696
Tool to read graphs from \e Nauty format data.
697 697
*/
698 698

	
699 699
/**
700 700
@defgroup concept Concepts
701 701
\brief Skeleton classes and concept checking classes
702 702

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

	
706 706
The purpose of the classes in this group is fourfold.
707 707

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

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

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

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

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

	
733 733
This group contains the skeletons and concept checking classes of
734 734
graph structures.
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2010
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 69
  /// \ref min_cost_flow "minimum cost flow" \ref amo93networkflows,
70 70
  /// \ref edmondskarp72theoretical. It is an efficient dual
71 71
  /// solution method.
72 72
  ///
73 73
  /// Most of the parameters of the problem (except for the digraph)
74 74
  /// can be given using separate functions, and the algorithm can be
75 75
  /// executed using the \ref run() function. If some parameters are not
76 76
  /// specified, then default values will be used.
77 77
  ///
78 78
  /// \tparam GR The digraph type the algorithm runs on.
79 79
  /// \tparam V The number type used for flow amounts, capacity bounds
80 80
  /// and supply values in the algorithm. By default, it is \c int.
81 81
  /// \tparam C The number type used for costs and potentials in the
82 82
  /// algorithm. By default, it is the same as \c V.
83 83
  /// \tparam TR The traits class that defines various types used by the
84 84
  /// algorithm. By default, it is \ref CapacityScalingDefaultTraits
85 85
  /// "CapacityScalingDefaultTraits<GR, V, C>".
86 86
  /// In most cases, this parameter should not be set directly,
87 87
  /// consider to use the named template parameters instead.
88 88
  ///
89 89
  /// \warning Both number types must be signed and all input data must
90 90
  /// be integer.
91
  /// \warning This algorithm does not support negative costs for such
92
  /// arcs that have infinite upper bound.
91
  /// \warning This algorithm does not support negative costs for
92
  /// arcs having infinite upper bound.
93 93
#ifdef DOXYGEN
94 94
  template <typename GR, typename V, typename C, typename TR>
95 95
#else
96 96
  template < typename GR, typename V = int, typename C = V,
97 97
             typename TR = CapacityScalingDefaultTraits<GR, V, C> >
98 98
#endif
99 99
  class CapacityScaling
100 100
  {
101 101
  public:
102 102

	
103 103
    /// The type of the digraph
104 104
    typedef typename TR::Digraph Digraph;
105 105
    /// The type of the flow amounts, capacity bounds and supply values
106 106
    typedef typename TR::Value Value;
107 107
    /// The type of the arc costs
108 108
    typedef typename TR::Cost Cost;
109 109

	
110 110
    /// The type of the heap used for internal Dijkstra computations
111 111
    typedef typename TR::Heap Heap;
112 112

	
113 113
    /// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm
114 114
    typedef TR Traits;
115 115

	
116 116
  public:
117 117

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

	
137 137
  private:
138 138

	
139 139
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
140 140

	
141 141
    typedef std::vector<int> IntVector;
142 142
    typedef std::vector<Value> ValueVector;
143 143
    typedef std::vector<Cost> CostVector;
144 144
    typedef std::vector<char> BoolVector;
145 145
    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
146 146

	
147 147
  private:
148 148

	
149 149
    // Data related to the underlying digraph
150 150
    const GR &_graph;
151 151
    int _node_num;
152 152
    int _arc_num;
153 153
    int _res_arc_num;
154 154
    int _root;
155 155

	
156 156
    // Parameters of the problem
157 157
    bool _have_lower;
158 158
    Value _sum_supply;
159 159

	
160 160
    // Data structures for storing the digraph
161 161
    IntNodeMap _node_id;
162 162
    IntArcMap _arc_idf;
163 163
    IntArcMap _arc_idb;
164 164
    IntVector _first_out;
165 165
    BoolVector _forward;
166 166
    IntVector _source;
167 167
    IntVector _target;
168 168
    IntVector _reverse;
169 169

	
170 170
    // Node and arc data
171 171
    ValueVector _lower;
172 172
    ValueVector _upper;
173 173
    CostVector _cost;
174 174
    ValueVector _supply;
175 175

	
176 176
    ValueVector _res_cap;
177 177
    CostVector _pi;
178 178
    ValueVector _excess;
179 179
    IntVector _excess_nodes;
180 180
    IntVector _deficit_nodes;
181 181

	
182 182
    Value _delta;
183 183
    int _factor;
184 184
    IntVector _pred;
185 185

	
186 186
  public:
187 187

	
188 188
    /// \brief Constant for infinite upper bounds (capacities).
189 189
    ///
190 190
    /// Constant for infinite upper bounds (capacities).
191 191
    /// It is \c std::numeric_limits<Value>::infinity() if available,
192 192
    /// \c std::numeric_limits<Value>::max() otherwise.
193 193
    const Value INF;
194 194

	
195 195
  private:
196 196

	
197 197
    // Special implementation of the Dijkstra algorithm for finding
198 198
    // shortest paths in the residual network of the digraph with
199 199
    // respect to the reduced arc costs and modifying the node
200 200
    // potentials according to the found distance labels.
201 201
    class ResidualDijkstra
202 202
    {
203 203
    private:
204 204

	
205 205
      int _node_num;
206 206
      bool _geq;
207 207
      const IntVector &_first_out;
208 208
      const IntVector &_target;
209 209
      const CostVector &_cost;
210 210
      const ValueVector &_res_cap;
211 211
      const ValueVector &_excess;
212 212
      CostVector &_pi;
213 213
      IntVector &_pred;
214 214

	
215 215
      IntVector _proc_nodes;
216 216
      CostVector _dist;
217 217

	
218 218
    public:
219 219

	
220 220
      ResidualDijkstra(CapacityScaling& cs) :
221 221
        _node_num(cs._node_num), _geq(cs._sum_supply < 0),
222 222
        _first_out(cs._first_out), _target(cs._target), _cost(cs._cost),
223 223
        _res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi),
224 224
        _pred(cs._pred), _dist(cs._node_num)
225 225
      {}
226 226

	
227 227
      int run(int s, Value delta = 1) {
228 228
        RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP);
229 229
        Heap heap(heap_cross_ref);
230 230
        heap.push(s, 0);
231 231
        _pred[s] = -1;
232 232
        _proc_nodes.clear();
233 233

	
234 234
        // Process nodes
235 235
        while (!heap.empty() && _excess[heap.top()] > -delta) {
236 236
          int u = heap.top(), v;
237 237
          Cost d = heap.prio() + _pi[u], dn;
238 238
          _dist[u] = heap.prio();
239 239
          _proc_nodes.push_back(u);
240 240
          heap.pop();
241 241

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

	
266 266
        // Update potentials of processed nodes
267 267
        int t = heap.top();
268 268
        Cost dt = heap.prio();
269 269
        for (int i = 0; i < int(_proc_nodes.size()); ++i) {
270 270
          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;
271 271
        }
272 272

	
273 273
        return t;
274 274
      }
275 275

	
276 276
    }; //class ResidualDijkstra
277 277

	
278 278
  public:
279 279

	
280 280
    /// \name Named Template Parameters
281 281
    /// @{
282 282

	
283 283
    template <typename T>
284 284
    struct SetHeapTraits : public Traits {
285 285
      typedef T Heap;
286 286
    };
287 287

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

	
302 302
    /// @}
303 303

	
304 304
  protected:
305 305

	
306 306
    CapacityScaling() {}
307 307

	
308 308
  public:
309 309

	
310 310
    /// \brief Constructor.
311 311
    ///
312 312
    /// The constructor of the class.
313 313
    ///
314 314
    /// \param graph The digraph the algorithm runs on.
315 315
    CapacityScaling(const GR& graph) :
316 316
      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
317 317
      INF(std::numeric_limits<Value>::has_infinity ?
318 318
          std::numeric_limits<Value>::infinity() :
319 319
          std::numeric_limits<Value>::max())
320 320
    {
321 321
      // Check the number types
322 322
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
323 323
        "The flow type of CapacityScaling must be signed");
324 324
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
325 325
        "The cost type of CapacityScaling must be signed");
326 326

	
327 327
      // Reset data structures
328 328
      reset();
329 329
    }
330 330

	
331 331
    /// \name Parameters
332 332
    /// The parameters of the algorithm can be specified using these
333 333
    /// functions.
334 334

	
335 335
    /// @{
336 336

	
337 337
    /// \brief Set the lower bounds on the arcs.
338 338
    ///
339 339
    /// This function sets the lower bounds on the arcs.
340 340
    /// If it is not used before calling \ref run(), the lower bounds
341 341
    /// will be set to zero on all arcs.
342 342
    ///
343 343
    /// \param map An arc map storing the lower bounds.
344 344
    /// Its \c Value type must be convertible to the \c Value type
345 345
    /// of the algorithm.
346 346
    ///
347 347
    /// \return <tt>(*this)</tt>
348 348
    template <typename LowerMap>
349 349
    CapacityScaling& lowerMap(const LowerMap& map) {
350 350
      _have_lower = true;
351 351
      for (ArcIt a(_graph); a != INVALID; ++a) {
352 352
        _lower[_arc_idf[a]] = map[a];
353 353
        _lower[_arc_idb[a]] = map[a];
354 354
      }
355 355
      return *this;
356 356
    }
357 357

	
358 358
    /// \brief Set the upper bounds (capacities) on the arcs.
359 359
    ///
360 360
    /// This function sets the upper bounds (capacities) on the arcs.
361 361
    /// If it is not used before calling \ref run(), the upper bounds
362 362
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
363 363
    /// unbounded from above).
364 364
    ///
365 365
    /// \param map An arc map storing the upper bounds.
366 366
    /// Its \c Value type must be convertible to the \c Value type
367 367
    /// of the algorithm.
368 368
    ///
369 369
    /// \return <tt>(*this)</tt>
370 370
    template<typename UpperMap>
371 371
    CapacityScaling& upperMap(const UpperMap& map) {
372 372
      for (ArcIt a(_graph); a != INVALID; ++a) {
373 373
        _upper[_arc_idf[a]] = map[a];
374 374
      }
375 375
      return *this;
376 376
    }
377 377

	
378 378
    /// \brief Set the costs of the arcs.
379 379
    ///
380 380
    /// This function sets the costs of the arcs.
381 381
    /// If it is not used before calling \ref run(), the costs
382 382
    /// will be set to \c 1 on all arcs.
383 383
    ///
384 384
    /// \param map An arc map storing the costs.
385 385
    /// Its \c Value type must be convertible to the \c Cost type
386 386
    /// of the algorithm.
387 387
    ///
388 388
    /// \return <tt>(*this)</tt>
389 389
    template<typename CostMap>
390 390
    CapacityScaling& costMap(const CostMap& map) {
391 391
      for (ArcIt a(_graph); a != INVALID; ++a) {
392 392
        _cost[_arc_idf[a]] =  map[a];
393 393
        _cost[_arc_idb[a]] = -map[a];
394 394
      }
395 395
      return *this;
396 396
    }
397 397

	
398 398
    /// \brief Set the supply values of the nodes.
399 399
    ///
400 400
    /// This function sets the supply values of the nodes.
401 401
    /// If neither this function nor \ref stSupply() is used before
402 402
    /// calling \ref run(), the supply of each node will be set to zero.
403 403
    ///
404 404
    /// \param map A node map storing the supply values.
405 405
    /// Its \c Value type must be convertible to the \c Value type
406 406
    /// of the algorithm.
407 407
    ///
408 408
    /// \return <tt>(*this)</tt>
409 409
    template<typename SupplyMap>
410 410
    CapacityScaling& supplyMap(const SupplyMap& map) {
411 411
      for (NodeIt n(_graph); n != INVALID; ++n) {
412 412
        _supply[_node_id[n]] = map[n];
413 413
      }
414 414
      return *this;
415 415
    }
416 416

	
417 417
    /// \brief Set single source and target nodes and a supply value.
418 418
    ///
419 419
    /// This function sets a single source node and a single target node
420 420
    /// and the required flow value.
421 421
    /// If neither this function nor \ref supplyMap() is used before
422 422
    /// calling \ref run(), the supply of each node will be set to zero.
423 423
    ///
424 424
    /// Using this function has the same effect as using \ref supplyMap()
425
    /// with such a map in which \c k is assigned to \c s, \c -k is
425
    /// with a map in which \c k is assigned to \c s, \c -k is
426 426
    /// assigned to \c t and all other nodes have zero supply value.
427 427
    ///
428 428
    /// \param s The source node.
429 429
    /// \param t The target node.
430 430
    /// \param k The required amount of flow from node \c s to node \c t
431 431
    /// (i.e. the supply of \c s and the demand of \c t).
432 432
    ///
433 433
    /// \return <tt>(*this)</tt>
434 434
    CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
435 435
      for (int i = 0; i != _node_num; ++i) {
436 436
        _supply[i] = 0;
437 437
      }
438 438
      _supply[_node_id[s]] =  k;
439 439
      _supply[_node_id[t]] = -k;
440 440
      return *this;
441 441
    }
442 442

	
443 443
    /// @}
444 444

	
445 445
    /// \name Execution control
446 446
    /// The algorithm can be executed using \ref run().
447 447

	
448 448
    /// @{
449 449

	
450 450
    /// \brief Run the algorithm.
451 451
    ///
452 452
    /// This function runs the algorithm.
453 453
    /// The paramters can be specified using functions \ref lowerMap(),
454 454
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
455 455
    /// For example,
456 456
    /// \code
457 457
    ///   CapacityScaling<ListDigraph> cs(graph);
458 458
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
459 459
    ///     .supplyMap(sup).run();
460 460
    /// \endcode
461 461
    ///
462 462
    /// This function can be called more than once. All the given parameters
463 463
    /// are kept for the next call, unless \ref resetParams() or \ref reset()
464 464
    /// is used, thus only the modified parameters have to be set again.
465 465
    /// If the underlying digraph was also modified after the construction
466 466
    /// of the class (or the last \ref reset() call), then the \ref reset()
467 467
    /// function must be called.
468 468
    ///
469 469
    /// \param factor The capacity scaling factor. It must be larger than
470 470
    /// one to use scaling. If it is less or equal to one, then scaling
471 471
    /// will be disabled.
472 472
    ///
473 473
    /// \return \c INFEASIBLE if no feasible flow exists,
474 474
    /// \n \c OPTIMAL if the problem has optimal solution
475 475
    /// (i.e. it is feasible and bounded), and the algorithm has found
476 476
    /// optimal flow and node potentials (primal and dual solutions),
477 477
    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
478 478
    /// and infinite upper bound. It means that the objective function
479 479
    /// is unbounded on that arc, however, note that it could actually be
480 480
    /// bounded over the feasible flows, but this algroithm cannot handle
481 481
    /// these cases.
482 482
    ///
483 483
    /// \see ProblemType
484 484
    /// \see resetParams(), reset()
485 485
    ProblemType run(int factor = 4) {
486 486
      _factor = factor;
487 487
      ProblemType pt = init();
488 488
      if (pt != OPTIMAL) return pt;
489 489
      return start();
490 490
    }
491 491

	
492 492
    /// \brief Reset all the parameters that have been given before.
493 493
    ///
494 494
    /// This function resets all the paramaters that have been given
495 495
    /// before using functions \ref lowerMap(), \ref upperMap(),
496 496
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
497 497
    ///
498 498
    /// It is useful for multiple \ref run() calls. Basically, all the given
499 499
    /// parameters are kept for the next \ref run() call, unless
500 500
    /// \ref resetParams() or \ref reset() is used.
501 501
    /// If the underlying digraph was also modified after the construction
502 502
    /// of the class or the last \ref reset() call, then the \ref reset()
503 503
    /// function must be used, otherwise \ref resetParams() is sufficient.
504 504
    ///
505 505
    /// For example,
506 506
    /// \code
507 507
    ///   CapacityScaling<ListDigraph> cs(graph);
508 508
    ///
509 509
    ///   // First run
510 510
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
511 511
    ///     .supplyMap(sup).run();
512 512
    ///
513 513
    ///   // Run again with modified cost map (resetParams() is not called,
514 514
    ///   // so only the cost map have to be set again)
515 515
    ///   cost[e] += 100;
516 516
    ///   cs.costMap(cost).run();
517 517
    ///
518 518
    ///   // Run again from scratch using resetParams()
519 519
    ///   // (the lower bounds will be set to zero on all arcs)
520 520
    ///   cs.resetParams();
521 521
    ///   cs.upperMap(capacity).costMap(cost)
522 522
    ///     .supplyMap(sup).run();
523 523
    /// \endcode
524 524
    ///
525 525
    /// \return <tt>(*this)</tt>
526 526
    ///
527 527
    /// \see reset(), run()
528 528
    CapacityScaling& resetParams() {
529 529
      for (int i = 0; i != _node_num; ++i) {
530 530
        _supply[i] = 0;
531 531
      }
532 532
      for (int j = 0; j != _res_arc_num; ++j) {
533 533
        _lower[j] = 0;
534 534
        _upper[j] = INF;
535 535
        _cost[j] = _forward[j] ? 1 : -1;
536 536
      }
537 537
      _have_lower = false;
538 538
      return *this;
539 539
    }
540 540

	
541 541
    /// \brief Reset the internal data structures and all the parameters
542 542
    /// that have been given before.
543 543
    ///
544 544
    /// This function resets the internal data structures and all the
545 545
    /// paramaters that have been given before using functions \ref lowerMap(),
546 546
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
547 547
    ///
548 548
    /// It is useful for multiple \ref run() calls. Basically, all the given
549 549
    /// parameters are kept for the next \ref run() call, unless
550 550
    /// \ref resetParams() or \ref reset() is used.
551 551
    /// If the underlying digraph was also modified after the construction
552 552
    /// of the class or the last \ref reset() call, then the \ref reset()
553 553
    /// function must be used, otherwise \ref resetParams() is sufficient.
554 554
    ///
555 555
    /// See \ref resetParams() for examples.
556 556
    ///
557 557
    /// \return <tt>(*this)</tt>
558 558
    ///
559 559
    /// \see resetParams(), run()
560 560
    CapacityScaling& reset() {
561 561
      // Resize vectors
562 562
      _node_num = countNodes(_graph);
563 563
      _arc_num = countArcs(_graph);
564 564
      _res_arc_num = 2 * (_arc_num + _node_num);
565 565
      _root = _node_num;
566 566
      ++_node_num;
567 567

	
568 568
      _first_out.resize(_node_num + 1);
569 569
      _forward.resize(_res_arc_num);
570 570
      _source.resize(_res_arc_num);
571 571
      _target.resize(_res_arc_num);
572 572
      _reverse.resize(_res_arc_num);
573 573

	
574 574
      _lower.resize(_res_arc_num);
575 575
      _upper.resize(_res_arc_num);
576 576
      _cost.resize(_res_arc_num);
577 577
      _supply.resize(_node_num);
578 578

	
579 579
      _res_cap.resize(_res_arc_num);
580 580
      _pi.resize(_node_num);
581 581
      _excess.resize(_node_num);
582 582
      _pred.resize(_node_num);
583 583

	
584 584
      // Copy the graph
585 585
      int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
586 586
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
587 587
        _node_id[n] = i;
588 588
      }
589 589
      i = 0;
590 590
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
591 591
        _first_out[i] = j;
592 592
        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
593 593
          _arc_idf[a] = j;
594 594
          _forward[j] = true;
595 595
          _source[j] = i;
596 596
          _target[j] = _node_id[_graph.runningNode(a)];
597 597
        }
598 598
        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
599 599
          _arc_idb[a] = j;
600 600
          _forward[j] = false;
601 601
          _source[j] = i;
602 602
          _target[j] = _node_id[_graph.runningNode(a)];
603 603
        }
604 604
        _forward[j] = false;
605 605
        _source[j] = i;
606 606
        _target[j] = _root;
607 607
        _reverse[j] = k;
608 608
        _forward[k] = true;
609 609
        _source[k] = _root;
610 610
        _target[k] = i;
611 611
        _reverse[k] = j;
612 612
        ++j; ++k;
613 613
      }
614 614
      _first_out[i] = j;
615 615
      _first_out[_node_num] = k;
616 616
      for (ArcIt a(_graph); a != INVALID; ++a) {
617 617
        int fi = _arc_idf[a];
Ignore white space 6 line context
... ...
@@ -258,385 +258,385 @@
258 258
  }
259 259

	
260 260
  /// \brief Function to count the edges in the graph.
261 261
  ///
262 262
  /// This function counts the edges in the graph.
263 263
  /// The complexity of the function is <em>O</em>(<em>m</em>), but for some
264 264
  /// graph structures it is specialized to run in <em>O</em>(1).
265 265
  ///
266 266
  /// \note If the graph contains a \c edgeNum() member function and a
267 267
  /// \c EdgeNumTag tag then this function calls directly the member
268 268
  /// function to query the cardinality of the edge set.
269 269
  template <typename Graph>
270 270
  inline int countEdges(const Graph& g) {
271 271
    return _core_bits::CountEdgesSelector<Graph>::count(g);
272 272

	
273 273
  }
274 274

	
275 275

	
276 276
  template <typename Graph, typename DegIt>
277 277
  inline int countNodeDegree(const Graph& _g, const typename Graph::Node& _n) {
278 278
    int num = 0;
279 279
    for (DegIt it(_g, _n); it != INVALID; ++it) {
280 280
      ++num;
281 281
    }
282 282
    return num;
283 283
  }
284 284

	
285 285
  /// \brief Function to count the number of the out-arcs from node \c n.
286 286
  ///
287 287
  /// This function counts the number of the out-arcs from node \c n
288 288
  /// in the graph \c g.
289 289
  template <typename Graph>
290 290
  inline int countOutArcs(const Graph& g,  const typename Graph::Node& n) {
291 291
    return countNodeDegree<Graph, typename Graph::OutArcIt>(g, n);
292 292
  }
293 293

	
294 294
  /// \brief Function to count the number of the in-arcs to node \c n.
295 295
  ///
296 296
  /// This function counts the number of the in-arcs to node \c n
297 297
  /// in the graph \c g.
298 298
  template <typename Graph>
299 299
  inline int countInArcs(const Graph& g,  const typename Graph::Node& n) {
300 300
    return countNodeDegree<Graph, typename Graph::InArcIt>(g, n);
301 301
  }
302 302

	
303 303
  /// \brief Function to count the number of the inc-edges to node \c n.
304 304
  ///
305 305
  /// This function counts the number of the inc-edges to node \c n
306 306
  /// in the undirected graph \c g.
307 307
  template <typename Graph>
308 308
  inline int countIncEdges(const Graph& g,  const typename Graph::Node& n) {
309 309
    return countNodeDegree<Graph, typename Graph::IncEdgeIt>(g, n);
310 310
  }
311 311

	
312 312
  namespace _core_bits {
313 313

	
314 314
    template <typename Digraph, typename Item, typename RefMap>
315 315
    class MapCopyBase {
316 316
    public:
317 317
      virtual void copy(const Digraph& from, const RefMap& refMap) = 0;
318 318

	
319 319
      virtual ~MapCopyBase() {}
320 320
    };
321 321

	
322 322
    template <typename Digraph, typename Item, typename RefMap,
323 323
              typename FromMap, typename ToMap>
324 324
    class MapCopy : public MapCopyBase<Digraph, Item, RefMap> {
325 325
    public:
326 326

	
327 327
      MapCopy(const FromMap& map, ToMap& tmap)
328 328
        : _map(map), _tmap(tmap) {}
329 329

	
330 330
      virtual void copy(const Digraph& digraph, const RefMap& refMap) {
331 331
        typedef typename ItemSetTraits<Digraph, Item>::ItemIt ItemIt;
332 332
        for (ItemIt it(digraph); it != INVALID; ++it) {
333 333
          _tmap.set(refMap[it], _map[it]);
334 334
        }
335 335
      }
336 336

	
337 337
    private:
338 338
      const FromMap& _map;
339 339
      ToMap& _tmap;
340 340
    };
341 341

	
342 342
    template <typename Digraph, typename Item, typename RefMap, typename It>
343 343
    class ItemCopy : public MapCopyBase<Digraph, Item, RefMap> {
344 344
    public:
345 345

	
346 346
      ItemCopy(const Item& item, It& it) : _item(item), _it(it) {}
347 347

	
348 348
      virtual void copy(const Digraph&, const RefMap& refMap) {
349 349
        _it = refMap[_item];
350 350
      }
351 351

	
352 352
    private:
353 353
      Item _item;
354 354
      It& _it;
355 355
    };
356 356

	
357 357
    template <typename Digraph, typename Item, typename RefMap, typename Ref>
358 358
    class RefCopy : public MapCopyBase<Digraph, Item, RefMap> {
359 359
    public:
360 360

	
361 361
      RefCopy(Ref& map) : _map(map) {}
362 362

	
363 363
      virtual void copy(const Digraph& digraph, const RefMap& refMap) {
364 364
        typedef typename ItemSetTraits<Digraph, Item>::ItemIt ItemIt;
365 365
        for (ItemIt it(digraph); it != INVALID; ++it) {
366 366
          _map.set(it, refMap[it]);
367 367
        }
368 368
      }
369 369

	
370 370
    private:
371 371
      Ref& _map;
372 372
    };
373 373

	
374 374
    template <typename Digraph, typename Item, typename RefMap,
375 375
              typename CrossRef>
376 376
    class CrossRefCopy : public MapCopyBase<Digraph, Item, RefMap> {
377 377
    public:
378 378

	
379 379
      CrossRefCopy(CrossRef& cmap) : _cmap(cmap) {}
380 380

	
381 381
      virtual void copy(const Digraph& digraph, const RefMap& refMap) {
382 382
        typedef typename ItemSetTraits<Digraph, Item>::ItemIt ItemIt;
383 383
        for (ItemIt it(digraph); it != INVALID; ++it) {
384 384
          _cmap.set(refMap[it], it);
385 385
        }
386 386
      }
387 387

	
388 388
    private:
389 389
      CrossRef& _cmap;
390 390
    };
391 391

	
392 392
    template <typename Digraph, typename Enable = void>
393 393
    struct DigraphCopySelector {
394 394
      template <typename From, typename NodeRefMap, typename ArcRefMap>
395 395
      static void copy(const From& from, Digraph &to,
396 396
                       NodeRefMap& nodeRefMap, ArcRefMap& arcRefMap) {
397 397
        to.clear();
398 398
        for (typename From::NodeIt it(from); it != INVALID; ++it) {
399 399
          nodeRefMap[it] = to.addNode();
400 400
        }
401 401
        for (typename From::ArcIt it(from); it != INVALID; ++it) {
402 402
          arcRefMap[it] = to.addArc(nodeRefMap[from.source(it)],
403 403
                                    nodeRefMap[from.target(it)]);
404 404
        }
405 405
      }
406 406
    };
407 407

	
408 408
    template <typename Digraph>
409 409
    struct DigraphCopySelector<
410 410
      Digraph,
411 411
      typename enable_if<typename Digraph::BuildTag, void>::type>
412 412
    {
413 413
      template <typename From, typename NodeRefMap, typename ArcRefMap>
414 414
      static void copy(const From& from, Digraph &to,
415 415
                       NodeRefMap& nodeRefMap, ArcRefMap& arcRefMap) {
416 416
        to.build(from, nodeRefMap, arcRefMap);
417 417
      }
418 418
    };
419 419

	
420 420
    template <typename Graph, typename Enable = void>
421 421
    struct GraphCopySelector {
422 422
      template <typename From, typename NodeRefMap, typename EdgeRefMap>
423 423
      static void copy(const From& from, Graph &to,
424 424
                       NodeRefMap& nodeRefMap, EdgeRefMap& edgeRefMap) {
425 425
        to.clear();
426 426
        for (typename From::NodeIt it(from); it != INVALID; ++it) {
427 427
          nodeRefMap[it] = to.addNode();
428 428
        }
429 429
        for (typename From::EdgeIt it(from); it != INVALID; ++it) {
430 430
          edgeRefMap[it] = to.addEdge(nodeRefMap[from.u(it)],
431 431
                                      nodeRefMap[from.v(it)]);
432 432
        }
433 433
      }
434 434
    };
435 435

	
436 436
    template <typename Graph>
437 437
    struct GraphCopySelector<
438 438
      Graph,
439 439
      typename enable_if<typename Graph::BuildTag, void>::type>
440 440
    {
441 441
      template <typename From, typename NodeRefMap, typename EdgeRefMap>
442 442
      static void copy(const From& from, Graph &to,
443 443
                       NodeRefMap& nodeRefMap, EdgeRefMap& edgeRefMap) {
444 444
        to.build(from, nodeRefMap, edgeRefMap);
445 445
      }
446 446
    };
447 447

	
448 448
  }
449 449

	
450
  /// Check whether a graph is undirected.
450
  /// \brief Check whether a graph is undirected.
451 451
  ///
452 452
  /// This function returns \c true if the given graph is undirected.
453 453
#ifdef DOXYGEN
454 454
  template <typename GR>
455 455
  bool undirected(const GR& g) { return false; }
456 456
#else
457 457
  template <typename GR>
458 458
  typename enable_if<UndirectedTagIndicator<GR>, bool>::type
459 459
  undirected(const GR&) {
460 460
    return true;
461 461
  }
462 462
  template <typename GR>
463 463
  typename disable_if<UndirectedTagIndicator<GR>, bool>::type
464 464
  undirected(const GR&) {
465 465
    return false;
466 466
  }
467 467
#endif
468 468

	
469 469
  /// \brief Class to copy a digraph.
470 470
  ///
471 471
  /// Class to copy a digraph to another digraph (duplicate a digraph). The
472 472
  /// simplest way of using it is through the \c digraphCopy() function.
473 473
  ///
474 474
  /// This class not only make a copy of a digraph, but it can create
475 475
  /// references and cross references between the nodes and arcs of
476 476
  /// the two digraphs, and it can copy maps to use with the newly created
477 477
  /// digraph.
478 478
  ///
479 479
  /// To make a copy from a digraph, first an instance of DigraphCopy
480 480
  /// should be created, then the data belongs to the digraph should
481 481
  /// assigned to copy. In the end, the \c run() member should be
482 482
  /// called.
483 483
  ///
484 484
  /// The next code copies a digraph with several data:
485 485
  ///\code
486 486
  ///  DigraphCopy<OrigGraph, NewGraph> cg(orig_graph, new_graph);
487 487
  ///  // Create references for the nodes
488 488
  ///  OrigGraph::NodeMap<NewGraph::Node> nr(orig_graph);
489 489
  ///  cg.nodeRef(nr);
490 490
  ///  // Create cross references (inverse) for the arcs
491 491
  ///  NewGraph::ArcMap<OrigGraph::Arc> acr(new_graph);
492 492
  ///  cg.arcCrossRef(acr);
493 493
  ///  // Copy an arc map
494 494
  ///  OrigGraph::ArcMap<double> oamap(orig_graph);
495 495
  ///  NewGraph::ArcMap<double> namap(new_graph);
496 496
  ///  cg.arcMap(oamap, namap);
497 497
  ///  // Copy a node
498 498
  ///  OrigGraph::Node on;
499 499
  ///  NewGraph::Node nn;
500 500
  ///  cg.node(on, nn);
501 501
  ///  // Execute copying
502 502
  ///  cg.run();
503 503
  ///\endcode
504 504
  template <typename From, typename To>
505 505
  class DigraphCopy {
506 506
  private:
507 507

	
508 508
    typedef typename From::Node Node;
509 509
    typedef typename From::NodeIt NodeIt;
510 510
    typedef typename From::Arc Arc;
511 511
    typedef typename From::ArcIt ArcIt;
512 512

	
513 513
    typedef typename To::Node TNode;
514 514
    typedef typename To::Arc TArc;
515 515

	
516 516
    typedef typename From::template NodeMap<TNode> NodeRefMap;
517 517
    typedef typename From::template ArcMap<TArc> ArcRefMap;
518 518

	
519 519
  public:
520 520

	
521 521
    /// \brief Constructor of DigraphCopy.
522 522
    ///
523 523
    /// Constructor of DigraphCopy for copying the content of the
524 524
    /// \c from digraph into the \c to digraph.
525 525
    DigraphCopy(const From& from, To& to)
526 526
      : _from(from), _to(to) {}
527 527

	
528 528
    /// \brief Destructor of DigraphCopy
529 529
    ///
530 530
    /// Destructor of DigraphCopy.
531 531
    ~DigraphCopy() {
532 532
      for (int i = 0; i < int(_node_maps.size()); ++i) {
533 533
        delete _node_maps[i];
534 534
      }
535 535
      for (int i = 0; i < int(_arc_maps.size()); ++i) {
536 536
        delete _arc_maps[i];
537 537
      }
538 538

	
539 539
    }
540 540

	
541 541
    /// \brief Copy the node references into the given map.
542 542
    ///
543 543
    /// This function copies the node references into the given map.
544 544
    /// The parameter should be a map, whose key type is the Node type of
545 545
    /// the source digraph, while the value type is the Node type of the
546 546
    /// destination digraph.
547 547
    template <typename NodeRef>
548 548
    DigraphCopy& nodeRef(NodeRef& map) {
549 549
      _node_maps.push_back(new _core_bits::RefCopy<From, Node,
550 550
                           NodeRefMap, NodeRef>(map));
551 551
      return *this;
552 552
    }
553 553

	
554 554
    /// \brief Copy the node cross references into the given map.
555 555
    ///
556 556
    /// This function copies the node cross references (reverse references)
557 557
    /// into the given map. The parameter should be a map, whose key type
558 558
    /// is the Node type of the destination digraph, while the value type is
559 559
    /// the Node type of the source digraph.
560 560
    template <typename NodeCrossRef>
561 561
    DigraphCopy& nodeCrossRef(NodeCrossRef& map) {
562 562
      _node_maps.push_back(new _core_bits::CrossRefCopy<From, Node,
563 563
                           NodeRefMap, NodeCrossRef>(map));
564 564
      return *this;
565 565
    }
566 566

	
567 567
    /// \brief Make a copy of the given node map.
568 568
    ///
569 569
    /// This function makes a copy of the given node map for the newly
570 570
    /// created digraph.
571 571
    /// The key type of the new map \c tmap should be the Node type of the
572 572
    /// destination digraph, and the key type of the original map \c map
573 573
    /// should be the Node type of the source digraph.
574 574
    template <typename FromMap, typename ToMap>
575 575
    DigraphCopy& nodeMap(const FromMap& map, ToMap& tmap) {
576 576
      _node_maps.push_back(new _core_bits::MapCopy<From, Node,
577 577
                           NodeRefMap, FromMap, ToMap>(map, tmap));
578 578
      return *this;
579 579
    }
580 580

	
581 581
    /// \brief Make a copy of the given node.
582 582
    ///
583 583
    /// This function makes a copy of the given node.
584 584
    DigraphCopy& node(const Node& node, TNode& tnode) {
585 585
      _node_maps.push_back(new _core_bits::ItemCopy<From, Node,
586 586
                           NodeRefMap, TNode>(node, tnode));
587 587
      return *this;
588 588
    }
589 589

	
590 590
    /// \brief Copy the arc references into the given map.
591 591
    ///
592 592
    /// This function copies the arc references into the given map.
593 593
    /// The parameter should be a map, whose key type is the Arc type of
594 594
    /// the source digraph, while the value type is the Arc type of the
595 595
    /// destination digraph.
596 596
    template <typename ArcRef>
597 597
    DigraphCopy& arcRef(ArcRef& map) {
598 598
      _arc_maps.push_back(new _core_bits::RefCopy<From, Arc,
599 599
                          ArcRefMap, ArcRef>(map));
600 600
      return *this;
601 601
    }
602 602

	
603 603
    /// \brief Copy the arc cross references into the given map.
604 604
    ///
605 605
    /// This function copies the arc cross references (reverse references)
606 606
    /// into the given map. The parameter should be a map, whose key type
607 607
    /// is the Arc type of the destination digraph, while the value type is
608 608
    /// the Arc type of the source digraph.
609 609
    template <typename ArcCrossRef>
610 610
    DigraphCopy& arcCrossRef(ArcCrossRef& map) {
611 611
      _arc_maps.push_back(new _core_bits::CrossRefCopy<From, Arc,
612 612
                          ArcRefMap, ArcCrossRef>(map));
613 613
      return *this;
614 614
    }
615 615

	
616 616
    /// \brief Make a copy of the given arc map.
617 617
    ///
618 618
    /// This function makes a copy of the given arc map for the newly
619 619
    /// created digraph.
620 620
    /// The key type of the new map \c tmap should be the Arc type of the
621 621
    /// destination digraph, and the key type of the original map \c map
622 622
    /// should be the Arc type of the source digraph.
623 623
    template <typename FromMap, typename ToMap>
624 624
    DigraphCopy& arcMap(const FromMap& map, ToMap& tmap) {
625 625
      _arc_maps.push_back(new _core_bits::MapCopy<From, Arc,
626 626
                          ArcRefMap, FromMap, ToMap>(map, tmap));
627 627
      return *this;
628 628
    }
629 629

	
630 630
    /// \brief Make a copy of the given arc.
631 631
    ///
632 632
    /// This function makes a copy of the given arc.
633 633
    DigraphCopy& arc(const Arc& arc, TArc& tarc) {
634 634
      _arc_maps.push_back(new _core_bits::ItemCopy<From, Arc,
635 635
                          ArcRefMap, TArc>(arc, tarc));
636 636
      return *this;
637 637
    }
638 638

	
639 639
    /// \brief Execute copying.
640 640
    ///
641 641
    /// This function executes the copying of the digraph along with the
642 642
    /// copying of the assigned data.
Ignore white space 384 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-2010
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 93
  /// push/augment and relabel operations for finding a \ref min_cost_flow
94 94
  /// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation,
95 95
  /// \ref goldberg97efficient, \ref bunnagel98efficient.
96 96
  /// It is a highly efficient primal-dual solution method, which
97 97
  /// can be viewed as the generalization of the \ref Preflow
98 98
  /// "preflow push-relabel" algorithm for the maximum flow problem.
99 99
  ///
100
  /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest
101
  /// implementations available in LEMON for this problem.
102
  ///
100 103
  /// Most of the parameters of the problem (except for the digraph)
101 104
  /// can be given using separate functions, and the algorithm can be
102 105
  /// executed using the \ref run() function. If some parameters are not
103 106
  /// specified, then default values will be used.
104 107
  ///
105 108
  /// \tparam GR The digraph type the algorithm runs on.
106 109
  /// \tparam V The number type used for flow amounts, capacity bounds
107 110
  /// and supply values in the algorithm. By default, it is \c int.
108 111
  /// \tparam C The number type used for costs and potentials in the
109 112
  /// algorithm. By default, it is the same as \c V.
110 113
  /// \tparam TR The traits class that defines various types used by the
111 114
  /// algorithm. By default, it is \ref CostScalingDefaultTraits
112 115
  /// "CostScalingDefaultTraits<GR, V, C>".
113 116
  /// In most cases, this parameter should not be set directly,
114 117
  /// consider to use the named template parameters instead.
115 118
  ///
116 119
  /// \warning Both number types must be signed and all input data must
117 120
  /// be integer.
118
  /// \warning This algorithm does not support negative costs for such
119
  /// arcs that have infinite upper bound.
121
  /// \warning This algorithm does not support negative costs for
122
  /// arcs having infinite upper bound.
120 123
  ///
121 124
  /// \note %CostScaling provides three different internal methods,
122 125
  /// from which the most efficient one is used by default.
123 126
  /// For more information, see \ref Method.
124 127
#ifdef DOXYGEN
125 128
  template <typename GR, typename V, typename C, typename TR>
126 129
#else
127 130
  template < typename GR, typename V = int, typename C = V,
128 131
             typename TR = CostScalingDefaultTraits<GR, V, C> >
129 132
#endif
130 133
  class CostScaling
131 134
  {
132 135
  public:
133 136

	
134 137
    /// The type of the digraph
135 138
    typedef typename TR::Digraph Digraph;
136 139
    /// The type of the flow amounts, capacity bounds and supply values
137 140
    typedef typename TR::Value Value;
138 141
    /// The type of the arc costs
139 142
    typedef typename TR::Cost Cost;
140 143

	
141 144
    /// \brief The large cost type
142 145
    ///
143 146
    /// The large cost type used for internal computations.
144 147
    /// By default, it is \c long \c long if the \c Cost type is integer,
145 148
    /// otherwise it is \c double.
146 149
    typedef typename TR::LargeCost LargeCost;
147 150

	
148 151
    /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
149 152
    typedef TR Traits;
150 153

	
151 154
  public:
152 155

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

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

	
199 202
  private:
200 203

	
201 204
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
202 205

	
203 206
    typedef std::vector<int> IntVector;
204 207
    typedef std::vector<Value> ValueVector;
205 208
    typedef std::vector<Cost> CostVector;
206 209
    typedef std::vector<LargeCost> LargeCostVector;
207 210
    typedef std::vector<char> BoolVector;
208 211
    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
209 212

	
210 213
  private:
211 214

	
212 215
    template <typename KT, typename VT>
213 216
    class StaticVectorMap {
214 217
    public:
215 218
      typedef KT Key;
216 219
      typedef VT Value;
217 220

	
218 221
      StaticVectorMap(std::vector<Value>& v) : _v(v) {}
219 222

	
220 223
      const Value& operator[](const Key& key) const {
221 224
        return _v[StaticDigraph::id(key)];
222 225
      }
223 226

	
224 227
      Value& operator[](const Key& key) {
225 228
        return _v[StaticDigraph::id(key)];
226 229
      }
227 230

	
228 231
      void set(const Key& key, const Value& val) {
229 232
        _v[StaticDigraph::id(key)] = val;
230 233
      }
231 234

	
232 235
    private:
233 236
      std::vector<Value>& _v;
234 237
    };
235 238

	
236 239
    typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
237 240
    typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
238 241

	
239 242
  private:
240 243

	
241 244
    // Data related to the underlying digraph
242 245
    const GR &_graph;
243 246
    int _node_num;
244 247
    int _arc_num;
245 248
    int _res_node_num;
246 249
    int _res_arc_num;
247 250
    int _root;
248 251

	
249 252
    // Parameters of the problem
250 253
    bool _have_lower;
251 254
    Value _sum_supply;
252 255
    int _sup_node_num;
253 256

	
254 257
    // Data structures for storing the digraph
255 258
    IntNodeMap _node_id;
256 259
    IntArcMap _arc_idf;
257 260
    IntArcMap _arc_idb;
258 261
    IntVector _first_out;
259 262
    BoolVector _forward;
260 263
    IntVector _source;
261 264
    IntVector _target;
262 265
    IntVector _reverse;
263 266

	
264 267
    // Node and arc data
265 268
    ValueVector _lower;
266 269
    ValueVector _upper;
267 270
    CostVector _scost;
268 271
    ValueVector _supply;
269 272

	
270 273
    ValueVector _res_cap;
271 274
    LargeCostVector _cost;
272 275
    LargeCostVector _pi;
273 276
    ValueVector _excess;
274 277
    IntVector _next_out;
275 278
    std::deque<int> _active_nodes;
276 279

	
277 280
    // Data for scaling
278 281
    LargeCost _epsilon;
279 282
    int _alpha;
280 283

	
281 284
    IntVector _buckets;
282 285
    IntVector _bucket_next;
283 286
    IntVector _bucket_prev;
284 287
    IntVector _rank;
285 288
    int _max_rank;
286 289

	
287 290
    // Data for a StaticDigraph structure
288 291
    typedef std::pair<int, int> IntPair;
289 292
    StaticDigraph _sgr;
290 293
    std::vector<IntPair> _arc_vec;
291 294
    std::vector<LargeCost> _cost_vec;
292 295
    LargeCostArcMap _cost_map;
293 296
    LargeCostNodeMap _pi_map;
294 297

	
295 298
  public:
296 299

	
297 300
    /// \brief Constant for infinite upper bounds (capacities).
298 301
    ///
299 302
    /// Constant for infinite upper bounds (capacities).
300 303
    /// It is \c std::numeric_limits<Value>::infinity() if available,
301 304
    /// \c std::numeric_limits<Value>::max() otherwise.
302 305
    const Value INF;
303 306

	
304 307
  public:
305 308

	
306 309
    /// \name Named Template Parameters
307 310
    /// @{
308 311

	
309 312
    template <typename T>
310 313
    struct SetLargeCostTraits : public Traits {
311 314
      typedef T LargeCost;
312 315
    };
313 316

	
314 317
    /// \brief \ref named-templ-param "Named parameter" for setting
315 318
    /// \c LargeCost type.
316 319
    ///
317 320
    /// \ref named-templ-param "Named parameter" for setting \c LargeCost
318 321
    /// type, which is used for internal computations in the algorithm.
319 322
    /// \c Cost must be convertible to \c LargeCost.
320 323
    template <typename T>
321 324
    struct SetLargeCost
322 325
      : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
323 326
      typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
324 327
    };
325 328

	
326 329
    /// @}
327 330

	
328 331
  protected:
329 332

	
330 333
    CostScaling() {}
331 334

	
332 335
  public:
333 336

	
334 337
    /// \brief Constructor.
335 338
    ///
336 339
    /// The constructor of the class.
337 340
    ///
338 341
    /// \param graph The digraph the algorithm runs on.
339 342
    CostScaling(const GR& graph) :
340 343
      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
341 344
      _cost_map(_cost_vec), _pi_map(_pi),
342 345
      INF(std::numeric_limits<Value>::has_infinity ?
343 346
          std::numeric_limits<Value>::infinity() :
344 347
          std::numeric_limits<Value>::max())
345 348
    {
346 349
      // Check the number types
347 350
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
348 351
        "The flow type of CostScaling must be signed");
349 352
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
350 353
        "The cost type of CostScaling must be signed");
351 354

	
352 355
      // Reset data structures
353 356
      reset();
354 357
    }
355 358

	
356 359
    /// \name Parameters
357 360
    /// The parameters of the algorithm can be specified using these
358 361
    /// functions.
359 362

	
360 363
    /// @{
361 364

	
362 365
    /// \brief Set the lower bounds on the arcs.
363 366
    ///
364 367
    /// This function sets the lower bounds on the arcs.
365 368
    /// If it is not used before calling \ref run(), the lower bounds
366 369
    /// will be set to zero on all arcs.
367 370
    ///
368 371
    /// \param map An arc map storing the lower bounds.
369 372
    /// Its \c Value type must be convertible to the \c Value type
370 373
    /// of the algorithm.
371 374
    ///
372 375
    /// \return <tt>(*this)</tt>
373 376
    template <typename LowerMap>
374 377
    CostScaling& lowerMap(const LowerMap& map) {
375 378
      _have_lower = true;
376 379
      for (ArcIt a(_graph); a != INVALID; ++a) {
377 380
        _lower[_arc_idf[a]] = map[a];
378 381
        _lower[_arc_idb[a]] = map[a];
379 382
      }
380 383
      return *this;
381 384
    }
382 385

	
383 386
    /// \brief Set the upper bounds (capacities) on the arcs.
384 387
    ///
385 388
    /// This function sets the upper bounds (capacities) on the arcs.
386 389
    /// If it is not used before calling \ref run(), the upper bounds
387 390
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
388 391
    /// unbounded from above).
389 392
    ///
390 393
    /// \param map An arc map storing the upper bounds.
391 394
    /// Its \c Value type must be convertible to the \c Value type
392 395
    /// of the algorithm.
393 396
    ///
394 397
    /// \return <tt>(*this)</tt>
395 398
    template<typename UpperMap>
396 399
    CostScaling& upperMap(const UpperMap& map) {
397 400
      for (ArcIt a(_graph); a != INVALID; ++a) {
398 401
        _upper[_arc_idf[a]] = map[a];
399 402
      }
400 403
      return *this;
401 404
    }
402 405

	
403 406
    /// \brief Set the costs of the arcs.
404 407
    ///
405 408
    /// This function sets the costs of the arcs.
406 409
    /// If it is not used before calling \ref run(), the costs
407 410
    /// will be set to \c 1 on all arcs.
408 411
    ///
409 412
    /// \param map An arc map storing the costs.
410 413
    /// Its \c Value type must be convertible to the \c Cost type
411 414
    /// of the algorithm.
412 415
    ///
413 416
    /// \return <tt>(*this)</tt>
414 417
    template<typename CostMap>
415 418
    CostScaling& costMap(const CostMap& map) {
416 419
      for (ArcIt a(_graph); a != INVALID; ++a) {
417 420
        _scost[_arc_idf[a]] =  map[a];
418 421
        _scost[_arc_idb[a]] = -map[a];
419 422
      }
420 423
      return *this;
421 424
    }
422 425

	
423 426
    /// \brief Set the supply values of the nodes.
424 427
    ///
425 428
    /// This function sets the supply values of the nodes.
426 429
    /// If neither this function nor \ref stSupply() is used before
427 430
    /// calling \ref run(), the supply of each node will be set to zero.
428 431
    ///
429 432
    /// \param map A node map storing the supply values.
430 433
    /// Its \c Value type must be convertible to the \c Value type
431 434
    /// of the algorithm.
432 435
    ///
433 436
    /// \return <tt>(*this)</tt>
434 437
    template<typename SupplyMap>
435 438
    CostScaling& supplyMap(const SupplyMap& map) {
436 439
      for (NodeIt n(_graph); n != INVALID; ++n) {
437 440
        _supply[_node_id[n]] = map[n];
438 441
      }
439 442
      return *this;
440 443
    }
441 444

	
442 445
    /// \brief Set single source and target nodes and a supply value.
443 446
    ///
444 447
    /// This function sets a single source node and a single target node
445 448
    /// and the required flow value.
446 449
    /// If neither this function nor \ref supplyMap() is used before
447 450
    /// calling \ref run(), the supply of each node will be set to zero.
448 451
    ///
449 452
    /// Using this function has the same effect as using \ref supplyMap()
450
    /// with such a map in which \c k is assigned to \c s, \c -k is
453
    /// with a map in which \c k is assigned to \c s, \c -k is
451 454
    /// assigned to \c t and all other nodes have zero supply value.
452 455
    ///
453 456
    /// \param s The source node.
454 457
    /// \param t The target node.
455 458
    /// \param k The required amount of flow from node \c s to node \c t
456 459
    /// (i.e. the supply of \c s and the demand of \c t).
457 460
    ///
458 461
    /// \return <tt>(*this)</tt>
459 462
    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
460 463
      for (int i = 0; i != _res_node_num; ++i) {
461 464
        _supply[i] = 0;
462 465
      }
463 466
      _supply[_node_id[s]] =  k;
464 467
      _supply[_node_id[t]] = -k;
465 468
      return *this;
466 469
    }
467 470

	
468 471
    /// @}
469 472

	
470 473
    /// \name Execution control
471 474
    /// The algorithm can be executed using \ref run().
472 475

	
473 476
    /// @{
474 477

	
475 478
    /// \brief Run the algorithm.
476 479
    ///
477 480
    /// This function runs the algorithm.
478 481
    /// The paramters can be specified using functions \ref lowerMap(),
479 482
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
480 483
    /// For example,
481 484
    /// \code
482 485
    ///   CostScaling<ListDigraph> cs(graph);
483 486
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
484 487
    ///     .supplyMap(sup).run();
485 488
    /// \endcode
486 489
    ///
487 490
    /// This function can be called more than once. All the given parameters
488 491
    /// are kept for the next call, unless \ref resetParams() or \ref reset()
489 492
    /// is used, thus only the modified parameters have to be set again.
490 493
    /// If the underlying digraph was also modified after the construction
491 494
    /// of the class (or the last \ref reset() call), then the \ref reset()
492 495
    /// function must be called.
493 496
    ///
494 497
    /// \param method The internal method that will be used in the
495 498
    /// algorithm. For more information, see \ref Method.
496 499
    /// \param factor The cost scaling factor. It must be larger than one.
497 500
    ///
498 501
    /// \return \c INFEASIBLE if no feasible flow exists,
499 502
    /// \n \c OPTIMAL if the problem has optimal solution
500 503
    /// (i.e. it is feasible and bounded), and the algorithm has found
501 504
    /// optimal flow and node potentials (primal and dual solutions),
502 505
    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
503 506
    /// and infinite upper bound. It means that the objective function
504 507
    /// is unbounded on that arc, however, note that it could actually be
505 508
    /// bounded over the feasible flows, but this algroithm cannot handle
506 509
    /// these cases.
507 510
    ///
508 511
    /// \see ProblemType, Method
509 512
    /// \see resetParams(), reset()
510 513
    ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
511 514
      _alpha = factor;
512 515
      ProblemType pt = init();
513 516
      if (pt != OPTIMAL) return pt;
514 517
      start(method);
515 518
      return OPTIMAL;
516 519
    }
517 520

	
518 521
    /// \brief Reset all the parameters that have been given before.
519 522
    ///
520 523
    /// This function resets all the paramaters that have been given
521 524
    /// before using functions \ref lowerMap(), \ref upperMap(),
522 525
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
523 526
    ///
524 527
    /// It is useful for multiple \ref run() calls. Basically, all the given
525 528
    /// parameters are kept for the next \ref run() call, unless
526 529
    /// \ref resetParams() or \ref reset() is used.
527 530
    /// If the underlying digraph was also modified after the construction
528 531
    /// of the class or the last \ref reset() call, then the \ref reset()
529 532
    /// function must be used, otherwise \ref resetParams() is sufficient.
530 533
    ///
531 534
    /// For example,
532 535
    /// \code
533 536
    ///   CostScaling<ListDigraph> cs(graph);
534 537
    ///
535 538
    ///   // First run
536 539
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
537 540
    ///     .supplyMap(sup).run();
538 541
    ///
539 542
    ///   // Run again with modified cost map (resetParams() is not called,
540 543
    ///   // so only the cost map have to be set again)
541 544
    ///   cost[e] += 100;
542 545
    ///   cs.costMap(cost).run();
543 546
    ///
544 547
    ///   // Run again from scratch using resetParams()
545 548
    ///   // (the lower bounds will be set to zero on all arcs)
546 549
    ///   cs.resetParams();
547 550
    ///   cs.upperMap(capacity).costMap(cost)
548 551
    ///     .supplyMap(sup).run();
549 552
    /// \endcode
550 553
    ///
551 554
    /// \return <tt>(*this)</tt>
552 555
    ///
553 556
    /// \see reset(), run()
554 557
    CostScaling& resetParams() {
555 558
      for (int i = 0; i != _res_node_num; ++i) {
556 559
        _supply[i] = 0;
557 560
      }
558 561
      int limit = _first_out[_root];
559 562
      for (int j = 0; j != limit; ++j) {
560 563
        _lower[j] = 0;
561 564
        _upper[j] = INF;
562 565
        _scost[j] = _forward[j] ? 1 : -1;
563 566
      }
564 567
      for (int j = limit; j != _res_arc_num; ++j) {
565 568
        _lower[j] = 0;
566 569
        _upper[j] = INF;
567 570
        _scost[j] = 0;
568 571
        _scost[_reverse[j]] = 0;
569 572
      }
570 573
      _have_lower = false;
571 574
      return *this;
572 575
    }
573 576

	
574 577
    /// \brief Reset all the parameters that have been given before.
575 578
    ///
576 579
    /// This function resets all the paramaters that have been given
577 580
    /// before using functions \ref lowerMap(), \ref upperMap(),
578 581
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
579 582
    ///
580 583
    /// It is useful for multiple run() calls. If this function is not
581 584
    /// used, all the parameters given before are kept for the next
582 585
    /// \ref run() call.
583 586
    /// However, the underlying digraph must not be modified after this
584 587
    /// class have been constructed, since it copies and extends the graph.
585 588
    /// \return <tt>(*this)</tt>
586 589
    CostScaling& reset() {
587 590
      // Resize vectors
588 591
      _node_num = countNodes(_graph);
589 592
      _arc_num = countArcs(_graph);
590 593
      _res_node_num = _node_num + 1;
591 594
      _res_arc_num = 2 * (_arc_num + _node_num);
592 595
      _root = _node_num;
593 596

	
594 597
      _first_out.resize(_res_node_num + 1);
595 598
      _forward.resize(_res_arc_num);
596 599
      _source.resize(_res_arc_num);
597 600
      _target.resize(_res_arc_num);
598 601
      _reverse.resize(_res_arc_num);
599 602

	
600 603
      _lower.resize(_res_arc_num);
601 604
      _upper.resize(_res_arc_num);
602 605
      _scost.resize(_res_arc_num);
603 606
      _supply.resize(_res_node_num);
604 607

	
605 608
      _res_cap.resize(_res_arc_num);
606 609
      _cost.resize(_res_arc_num);
607 610
      _pi.resize(_res_node_num);
608 611
      _excess.resize(_res_node_num);
609 612
      _next_out.resize(_res_node_num);
610 613

	
611 614
      _arc_vec.reserve(_res_arc_num);
612 615
      _cost_vec.reserve(_res_arc_num);
613 616

	
614 617
      // Copy the graph
615 618
      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
616 619
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
617 620
        _node_id[n] = i;
618 621
      }
619 622
      i = 0;
620 623
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
621 624
        _first_out[i] = j;
622 625
        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
623 626
          _arc_idf[a] = j;
624 627
          _forward[j] = true;
625 628
          _source[j] = i;
626 629
          _target[j] = _node_id[_graph.runningNode(a)];
627 630
        }
628 631
        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
629 632
          _arc_idb[a] = j;
630 633
          _forward[j] = false;
631 634
          _source[j] = i;
632 635
          _target[j] = _node_id[_graph.runningNode(a)];
633 636
        }
634 637
        _forward[j] = false;
635 638
        _source[j] = i;
636 639
        _target[j] = _root;
637 640
        _reverse[j] = k;
638 641
        _forward[k] = true;
639 642
        _source[k] = _root;
640 643
        _target[k] = i;
641 644
        _reverse[k] = j;
642 645
        ++j; ++k;
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2010
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_CYCLE_CANCELING_H
20 20
#define LEMON_CYCLE_CANCELING_H
21 21

	
22 22
/// \ingroup min_cost_flow_algs
23 23
/// \file
24 24
/// \brief Cycle-canceling algorithms for finding a minimum cost flow.
25 25

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

	
29 29
#include <lemon/core.h>
30 30
#include <lemon/maps.h>
31 31
#include <lemon/path.h>
32 32
#include <lemon/math.h>
33 33
#include <lemon/static_graph.h>
34 34
#include <lemon/adaptors.h>
35 35
#include <lemon/circulation.h>
36 36
#include <lemon/bellman_ford.h>
37 37
#include <lemon/howard_mmc.h>
38 38

	
39 39
namespace lemon {
40 40

	
41 41
  /// \addtogroup min_cost_flow_algs
42 42
  /// @{
43 43

	
44 44
  /// \brief Implementation of cycle-canceling algorithms for
45 45
  /// finding a \ref min_cost_flow "minimum cost flow".
46 46
  ///
47 47
  /// \ref CycleCanceling implements three different cycle-canceling
48 48
  /// algorithms for finding a \ref min_cost_flow "minimum cost flow"
49 49
  /// \ref amo93networkflows, \ref klein67primal,
50 50
  /// \ref goldberg89cyclecanceling.
51 51
  /// The most efficent one (both theoretically and practically)
52 52
  /// is the \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" algorithm,
53 53
  /// thus it is the default method.
54 54
  /// It is strongly polynomial, but in practice, it is typically much
55 55
  /// slower than the scaling algorithms and NetworkSimplex.
56 56
  ///
57 57
  /// Most of the parameters of the problem (except for the digraph)
58 58
  /// can be given using separate functions, and the algorithm can be
59 59
  /// executed using the \ref run() function. If some parameters are not
60 60
  /// specified, then default values will be used.
61 61
  ///
62 62
  /// \tparam GR The digraph type the algorithm runs on.
63 63
  /// \tparam V The number type used for flow amounts, capacity bounds
64 64
  /// and supply values in the algorithm. By default, it is \c int.
65 65
  /// \tparam C The number type used for costs and potentials in the
66 66
  /// algorithm. By default, it is the same as \c V.
67 67
  ///
68 68
  /// \warning Both number types must be signed and all input data must
69 69
  /// be integer.
70
  /// \warning This algorithm does not support negative costs for such
71
  /// arcs that have infinite upper bound.
70
  /// \warning This algorithm does not support negative costs for
71
  /// arcs having infinite upper bound.
72 72
  ///
73 73
  /// \note For more information about the three available methods,
74 74
  /// see \ref Method.
75 75
#ifdef DOXYGEN
76 76
  template <typename GR, typename V, typename C>
77 77
#else
78 78
  template <typename GR, typename V = int, typename C = V>
79 79
#endif
80 80
  class CycleCanceling
81 81
  {
82 82
  public:
83 83

	
84 84
    /// The type of the digraph
85 85
    typedef GR Digraph;
86 86
    /// The type of the flow amounts, capacity bounds and supply values
87 87
    typedef V Value;
88 88
    /// The type of the arc costs
89 89
    typedef C Cost;
90 90

	
91 91
  public:
92 92

	
93 93
    /// \brief Problem type constants for the \c run() function.
94 94
    ///
95 95
    /// Enum type containing the problem type constants that can be
96 96
    /// returned by the \ref run() function of the algorithm.
97 97
    enum ProblemType {
98 98
      /// The problem has no feasible solution (flow).
99 99
      INFEASIBLE,
100 100
      /// The problem has optimal solution (i.e. it is feasible and
101 101
      /// bounded), and the algorithm has found optimal flow and node
102 102
      /// potentials (primal and dual solutions).
103 103
      OPTIMAL,
104 104
      /// The digraph contains an arc of negative cost and infinite
105 105
      /// upper bound. It means that the objective function is unbounded
106 106
      /// on that arc, however, note that it could actually be bounded
107 107
      /// over the feasible flows, but this algroithm cannot handle
108 108
      /// these cases.
109 109
      UNBOUNDED
110 110
    };
111 111

	
112 112
    /// \brief Constants for selecting the used method.
113 113
    ///
114 114
    /// Enum type containing constants for selecting the used method
115 115
    /// for the \ref run() function.
116 116
    ///
117 117
    /// \ref CycleCanceling provides three different cycle-canceling
118 118
    /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten"
119
    /// is used, which proved to be the most efficient and the most robust
120
    /// on various test inputs.
119
    /// is used, which is by far the most efficient and the most robust.
121 120
    /// However, the other methods can be selected using the \ref run()
122 121
    /// function with the proper parameter.
123 122
    enum Method {
124 123
      /// A simple cycle-canceling method, which uses the
125 124
      /// \ref BellmanFord "Bellman-Ford" algorithm with limited iteration
126 125
      /// number for detecting negative cycles in the residual network.
127 126
      SIMPLE_CYCLE_CANCELING,
128 127
      /// The "Minimum Mean Cycle-Canceling" algorithm, which is a
129 128
      /// well-known strongly polynomial method
130 129
      /// \ref goldberg89cyclecanceling. It improves along a
131 130
      /// \ref min_mean_cycle "minimum mean cycle" in each iteration.
132 131
      /// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)).
133 132
      MINIMUM_MEAN_CYCLE_CANCELING,
134 133
      /// The "Cancel And Tighten" algorithm, which can be viewed as an
135 134
      /// improved version of the previous method
136 135
      /// \ref goldberg89cyclecanceling.
137 136
      /// It is faster both in theory and in practice, its running time
138 137
      /// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)).
139 138
      CANCEL_AND_TIGHTEN
140 139
    };
141 140

	
142 141
  private:
143 142

	
144 143
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
145 144

	
146 145
    typedef std::vector<int> IntVector;
147 146
    typedef std::vector<double> DoubleVector;
148 147
    typedef std::vector<Value> ValueVector;
149 148
    typedef std::vector<Cost> CostVector;
150 149
    typedef std::vector<char> BoolVector;
151 150
    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
152 151

	
153 152
  private:
154 153

	
155 154
    template <typename KT, typename VT>
156 155
    class StaticVectorMap {
157 156
    public:
158 157
      typedef KT Key;
159 158
      typedef VT Value;
160 159

	
161 160
      StaticVectorMap(std::vector<Value>& v) : _v(v) {}
162 161

	
163 162
      const Value& operator[](const Key& key) const {
164 163
        return _v[StaticDigraph::id(key)];
165 164
      }
166 165

	
167 166
      Value& operator[](const Key& key) {
168 167
        return _v[StaticDigraph::id(key)];
169 168
      }
170 169

	
171 170
      void set(const Key& key, const Value& val) {
172 171
        _v[StaticDigraph::id(key)] = val;
173 172
      }
174 173

	
175 174
    private:
176 175
      std::vector<Value>& _v;
177 176
    };
178 177

	
179 178
    typedef StaticVectorMap<StaticDigraph::Node, Cost> CostNodeMap;
180 179
    typedef StaticVectorMap<StaticDigraph::Arc, Cost> CostArcMap;
181 180

	
182 181
  private:
183 182

	
184 183

	
185 184
    // Data related to the underlying digraph
186 185
    const GR &_graph;
187 186
    int _node_num;
188 187
    int _arc_num;
189 188
    int _res_node_num;
190 189
    int _res_arc_num;
191 190
    int _root;
192 191

	
193 192
    // Parameters of the problem
194 193
    bool _have_lower;
195 194
    Value _sum_supply;
196 195

	
197 196
    // Data structures for storing the digraph
198 197
    IntNodeMap _node_id;
199 198
    IntArcMap _arc_idf;
200 199
    IntArcMap _arc_idb;
201 200
    IntVector _first_out;
202 201
    BoolVector _forward;
203 202
    IntVector _source;
204 203
    IntVector _target;
205 204
    IntVector _reverse;
206 205

	
207 206
    // Node and arc data
208 207
    ValueVector _lower;
209 208
    ValueVector _upper;
210 209
    CostVector _cost;
211 210
    ValueVector _supply;
212 211

	
213 212
    ValueVector _res_cap;
214 213
    CostVector _pi;
215 214

	
216 215
    // Data for a StaticDigraph structure
217 216
    typedef std::pair<int, int> IntPair;
218 217
    StaticDigraph _sgr;
219 218
    std::vector<IntPair> _arc_vec;
220 219
    std::vector<Cost> _cost_vec;
221 220
    IntVector _id_vec;
222 221
    CostArcMap _cost_map;
223 222
    CostNodeMap _pi_map;
224 223

	
225 224
  public:
226 225

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

	
234 233
  public:
235 234

	
236 235
    /// \brief Constructor.
237 236
    ///
238 237
    /// The constructor of the class.
239 238
    ///
240 239
    /// \param graph The digraph the algorithm runs on.
241 240
    CycleCanceling(const GR& graph) :
242 241
      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
243 242
      _cost_map(_cost_vec), _pi_map(_pi),
244 243
      INF(std::numeric_limits<Value>::has_infinity ?
245 244
          std::numeric_limits<Value>::infinity() :
246 245
          std::numeric_limits<Value>::max())
247 246
    {
248 247
      // Check the number types
249 248
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
250 249
        "The flow type of CycleCanceling must be signed");
251 250
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
252 251
        "The cost type of CycleCanceling must be signed");
253 252

	
254 253
      // Reset data structures
255 254
      reset();
256 255
    }
257 256

	
258 257
    /// \name Parameters
259 258
    /// The parameters of the algorithm can be specified using these
260 259
    /// functions.
261 260

	
262 261
    /// @{
263 262

	
264 263
    /// \brief Set the lower bounds on the arcs.
265 264
    ///
266 265
    /// This function sets the lower bounds on the arcs.
267 266
    /// If it is not used before calling \ref run(), the lower bounds
268 267
    /// will be set to zero on all arcs.
269 268
    ///
270 269
    /// \param map An arc map storing the lower bounds.
271 270
    /// Its \c Value type must be convertible to the \c Value type
272 271
    /// of the algorithm.
273 272
    ///
274 273
    /// \return <tt>(*this)</tt>
275 274
    template <typename LowerMap>
276 275
    CycleCanceling& lowerMap(const LowerMap& map) {
277 276
      _have_lower = true;
278 277
      for (ArcIt a(_graph); a != INVALID; ++a) {
279 278
        _lower[_arc_idf[a]] = map[a];
280 279
        _lower[_arc_idb[a]] = map[a];
281 280
      }
282 281
      return *this;
283 282
    }
284 283

	
285 284
    /// \brief Set the upper bounds (capacities) on the arcs.
286 285
    ///
287 286
    /// This function sets the upper bounds (capacities) on the arcs.
288 287
    /// If it is not used before calling \ref run(), the upper bounds
289 288
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
290 289
    /// unbounded from above).
291 290
    ///
292 291
    /// \param map An arc map storing the upper bounds.
293 292
    /// Its \c Value type must be convertible to the \c Value type
294 293
    /// of the algorithm.
295 294
    ///
296 295
    /// \return <tt>(*this)</tt>
297 296
    template<typename UpperMap>
298 297
    CycleCanceling& upperMap(const UpperMap& map) {
299 298
      for (ArcIt a(_graph); a != INVALID; ++a) {
300 299
        _upper[_arc_idf[a]] = map[a];
301 300
      }
302 301
      return *this;
303 302
    }
304 303

	
305 304
    /// \brief Set the costs of the arcs.
306 305
    ///
307 306
    /// This function sets the costs of the arcs.
308 307
    /// If it is not used before calling \ref run(), the costs
309 308
    /// will be set to \c 1 on all arcs.
310 309
    ///
311 310
    /// \param map An arc map storing the costs.
312 311
    /// Its \c Value type must be convertible to the \c Cost type
313 312
    /// of the algorithm.
314 313
    ///
315 314
    /// \return <tt>(*this)</tt>
316 315
    template<typename CostMap>
317 316
    CycleCanceling& costMap(const CostMap& map) {
318 317
      for (ArcIt a(_graph); a != INVALID; ++a) {
319 318
        _cost[_arc_idf[a]] =  map[a];
320 319
        _cost[_arc_idb[a]] = -map[a];
321 320
      }
322 321
      return *this;
323 322
    }
324 323

	
325 324
    /// \brief Set the supply values of the nodes.
326 325
    ///
327 326
    /// This function sets the supply values of the nodes.
328 327
    /// If neither this function nor \ref stSupply() is used before
329 328
    /// calling \ref run(), the supply of each node will be set to zero.
330 329
    ///
331 330
    /// \param map A node map storing the supply values.
332 331
    /// Its \c Value type must be convertible to the \c Value type
333 332
    /// of the algorithm.
334 333
    ///
335 334
    /// \return <tt>(*this)</tt>
336 335
    template<typename SupplyMap>
337 336
    CycleCanceling& supplyMap(const SupplyMap& map) {
338 337
      for (NodeIt n(_graph); n != INVALID; ++n) {
339 338
        _supply[_node_id[n]] = map[n];
340 339
      }
341 340
      return *this;
342 341
    }
343 342

	
344 343
    /// \brief Set single source and target nodes and a supply value.
345 344
    ///
346 345
    /// This function sets a single source node and a single target node
347 346
    /// and the required flow value.
348 347
    /// If neither this function nor \ref supplyMap() is used before
349 348
    /// calling \ref run(), the supply of each node will be set to zero.
350 349
    ///
351 350
    /// Using this function has the same effect as using \ref supplyMap()
352
    /// with such a map in which \c k is assigned to \c s, \c -k is
351
    /// with a map in which \c k is assigned to \c s, \c -k is
353 352
    /// assigned to \c t and all other nodes have zero supply value.
354 353
    ///
355 354
    /// \param s The source node.
356 355
    /// \param t The target node.
357 356
    /// \param k The required amount of flow from node \c s to node \c t
358 357
    /// (i.e. the supply of \c s and the demand of \c t).
359 358
    ///
360 359
    /// \return <tt>(*this)</tt>
361 360
    CycleCanceling& stSupply(const Node& s, const Node& t, Value k) {
362 361
      for (int i = 0; i != _res_node_num; ++i) {
363 362
        _supply[i] = 0;
364 363
      }
365 364
      _supply[_node_id[s]] =  k;
366 365
      _supply[_node_id[t]] = -k;
367 366
      return *this;
368 367
    }
369 368

	
370 369
    /// @}
371 370

	
372 371
    /// \name Execution control
373 372
    /// The algorithm can be executed using \ref run().
374 373

	
375 374
    /// @{
376 375

	
377 376
    /// \brief Run the algorithm.
378 377
    ///
379 378
    /// This function runs the algorithm.
380 379
    /// The paramters can be specified using functions \ref lowerMap(),
381 380
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
382 381
    /// For example,
383 382
    /// \code
384 383
    ///   CycleCanceling<ListDigraph> cc(graph);
385 384
    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
386 385
    ///     .supplyMap(sup).run();
387 386
    /// \endcode
388 387
    ///
389 388
    /// This function can be called more than once. All the given parameters
390 389
    /// are kept for the next call, unless \ref resetParams() or \ref reset()
391 390
    /// is used, thus only the modified parameters have to be set again.
392 391
    /// If the underlying digraph was also modified after the construction
393 392
    /// of the class (or the last \ref reset() call), then the \ref reset()
394 393
    /// function must be called.
395 394
    ///
396 395
    /// \param method The cycle-canceling method that will be used.
397 396
    /// For more information, see \ref Method.
398 397
    ///
399 398
    /// \return \c INFEASIBLE if no feasible flow exists,
400 399
    /// \n \c OPTIMAL if the problem has optimal solution
401 400
    /// (i.e. it is feasible and bounded), and the algorithm has found
402 401
    /// optimal flow and node potentials (primal and dual solutions),
403 402
    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
404 403
    /// and infinite upper bound. It means that the objective function
405 404
    /// is unbounded on that arc, however, note that it could actually be
406 405
    /// bounded over the feasible flows, but this algroithm cannot handle
407 406
    /// these cases.
408 407
    ///
409 408
    /// \see ProblemType, Method
410 409
    /// \see resetParams(), reset()
411 410
    ProblemType run(Method method = CANCEL_AND_TIGHTEN) {
412 411
      ProblemType pt = init();
413 412
      if (pt != OPTIMAL) return pt;
414 413
      start(method);
415 414
      return OPTIMAL;
416 415
    }
417 416

	
418 417
    /// \brief Reset all the parameters that have been given before.
419 418
    ///
420 419
    /// This function resets all the paramaters that have been given
421 420
    /// before using functions \ref lowerMap(), \ref upperMap(),
422 421
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
423 422
    ///
424 423
    /// It is useful for multiple \ref run() calls. Basically, all the given
425 424
    /// parameters are kept for the next \ref run() call, unless
426 425
    /// \ref resetParams() or \ref reset() is used.
427 426
    /// If the underlying digraph was also modified after the construction
428 427
    /// of the class or the last \ref reset() call, then the \ref reset()
429 428
    /// function must be used, otherwise \ref resetParams() is sufficient.
430 429
    ///
431 430
    /// For example,
432 431
    /// \code
433 432
    ///   CycleCanceling<ListDigraph> cs(graph);
434 433
    ///
435 434
    ///   // First run
436 435
    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
437 436
    ///     .supplyMap(sup).run();
438 437
    ///
439 438
    ///   // Run again with modified cost map (resetParams() is not called,
440 439
    ///   // so only the cost map have to be set again)
441 440
    ///   cost[e] += 100;
442 441
    ///   cc.costMap(cost).run();
443 442
    ///
444 443
    ///   // Run again from scratch using resetParams()
445 444
    ///   // (the lower bounds will be set to zero on all arcs)
446 445
    ///   cc.resetParams();
447 446
    ///   cc.upperMap(capacity).costMap(cost)
448 447
    ///     .supplyMap(sup).run();
449 448
    /// \endcode
450 449
    ///
451 450
    /// \return <tt>(*this)</tt>
452 451
    ///
453 452
    /// \see reset(), run()
454 453
    CycleCanceling& resetParams() {
455 454
      for (int i = 0; i != _res_node_num; ++i) {
456 455
        _supply[i] = 0;
457 456
      }
458 457
      int limit = _first_out[_root];
459 458
      for (int j = 0; j != limit; ++j) {
460 459
        _lower[j] = 0;
461 460
        _upper[j] = INF;
462 461
        _cost[j] = _forward[j] ? 1 : -1;
463 462
      }
464 463
      for (int j = limit; j != _res_arc_num; ++j) {
465 464
        _lower[j] = 0;
466 465
        _upper[j] = INF;
467 466
        _cost[j] = 0;
468 467
        _cost[_reverse[j]] = 0;
469 468
      }
470 469
      _have_lower = false;
471 470
      return *this;
472 471
    }
473 472

	
474 473
    /// \brief Reset the internal data structures and all the parameters
475 474
    /// that have been given before.
476 475
    ///
477 476
    /// This function resets the internal data structures and all the
478 477
    /// paramaters that have been given before using functions \ref lowerMap(),
479 478
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
480 479
    ///
481 480
    /// It is useful for multiple \ref run() calls. Basically, all the given
482 481
    /// parameters are kept for the next \ref run() call, unless
483 482
    /// \ref resetParams() or \ref reset() is used.
484 483
    /// If the underlying digraph was also modified after the construction
485 484
    /// of the class or the last \ref reset() call, then the \ref reset()
486 485
    /// function must be used, otherwise \ref resetParams() is sufficient.
487 486
    ///
488 487
    /// See \ref resetParams() for examples.
489 488
    ///
490 489
    /// \return <tt>(*this)</tt>
491 490
    ///
492 491
    /// \see resetParams(), run()
493 492
    CycleCanceling& reset() {
494 493
      // Resize vectors
495 494
      _node_num = countNodes(_graph);
496 495
      _arc_num = countArcs(_graph);
497 496
      _res_node_num = _node_num + 1;
498 497
      _res_arc_num = 2 * (_arc_num + _node_num);
499 498
      _root = _node_num;
500 499

	
501 500
      _first_out.resize(_res_node_num + 1);
502 501
      _forward.resize(_res_arc_num);
503 502
      _source.resize(_res_arc_num);
504 503
      _target.resize(_res_arc_num);
505 504
      _reverse.resize(_res_arc_num);
506 505

	
507 506
      _lower.resize(_res_arc_num);
508 507
      _upper.resize(_res_arc_num);
509 508
      _cost.resize(_res_arc_num);
510 509
      _supply.resize(_res_node_num);
511 510

	
512 511
      _res_cap.resize(_res_arc_num);
513 512
      _pi.resize(_res_node_num);
514 513

	
515 514
      _arc_vec.reserve(_res_arc_num);
516 515
      _cost_vec.reserve(_res_arc_num);
517 516
      _id_vec.reserve(_res_arc_num);
518 517

	
519 518
      // Copy the graph
520 519
      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
521 520
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
522 521
        _node_id[n] = i;
523 522
      }
524 523
      i = 0;
525 524
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
526 525
        _first_out[i] = j;
527 526
        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
528 527
          _arc_idf[a] = j;
529 528
          _forward[j] = true;
530 529
          _source[j] = i;
531 530
          _target[j] = _node_id[_graph.runningNode(a)];
532 531
        }
533 532
        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
534 533
          _arc_idb[a] = j;
535 534
          _forward[j] = false;
536 535
          _source[j] = i;
537 536
          _target[j] = _node_id[_graph.runningNode(a)];
538 537
        }
539 538
        _forward[j] = false;
540 539
        _source[j] = i;
541 540
        _target[j] = _root;
542 541
        _reverse[j] = k;
543 542
        _forward[k] = true;
544 543
        _source[k] = _root;
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2010
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_EULER_H
20 20
#define LEMON_EULER_H
21 21

	
22 22
#include<lemon/core.h>
23 23
#include<lemon/adaptors.h>
24 24
#include<lemon/connectivity.h>
25 25
#include <list>
26 26

	
27 27
/// \ingroup graph_properties
28 28
/// \file
29 29
/// \brief Euler tour iterators and a function for checking the \e Eulerian
30 30
/// property.
31 31
///
32 32
///This file provides Euler tour iterators and a function to check
33 33
///if a (di)graph is \e Eulerian.
34 34

	
35 35
namespace lemon {
36 36

	
37 37
  ///Euler tour iterator for digraphs.
38 38

	
39
  /// \ingroup graph_prop
39
  /// \ingroup graph_properties
40 40
  ///This iterator provides an Euler tour (Eulerian circuit) of a \e directed
41 41
  ///graph (if there exists) and it converts to the \c Arc type of the digraph.
42 42
  ///
43 43
  ///For example, if the given digraph has an Euler tour (i.e it has only one
44 44
  ///non-trivial component and the in-degree is equal to the out-degree
45 45
  ///for all nodes), then the following code will put the arcs of \c g
46 46
  ///to the vector \c et according to an Euler tour of \c g.
47 47
  ///\code
48 48
  ///  std::vector<ListDigraph::Arc> et;
49 49
  ///  for(DiEulerIt<ListDigraph> e(g); e!=INVALID; ++e)
50 50
  ///    et.push_back(e);
51 51
  ///\endcode
52 52
  ///If \c g has no Euler tour, then the resulted walk will not be closed
53 53
  ///or not contain all arcs.
54 54
  ///\sa EulerIt
55 55
  template<typename GR>
56 56
  class DiEulerIt
57 57
  {
58 58
    typedef typename GR::Node Node;
59 59
    typedef typename GR::NodeIt NodeIt;
60 60
    typedef typename GR::Arc Arc;
61 61
    typedef typename GR::ArcIt ArcIt;
62 62
    typedef typename GR::OutArcIt OutArcIt;
63 63
    typedef typename GR::InArcIt InArcIt;
64 64

	
65 65
    const GR &g;
66 66
    typename GR::template NodeMap<OutArcIt> narc;
67 67
    std::list<Arc> euler;
68 68

	
69 69
  public:
70 70

	
71 71
    ///Constructor
72 72

	
73 73
    ///Constructor.
74 74
    ///\param gr A digraph.
75 75
    ///\param start The starting point of the tour. If it is not given,
76 76
    ///the tour will start from the first node that has an outgoing arc.
77 77
    DiEulerIt(const GR &gr, typename GR::Node start = INVALID)
78 78
      : g(gr), narc(g)
79 79
    {
80 80
      if (start==INVALID) {
81 81
        NodeIt n(g);
82 82
        while (n!=INVALID && OutArcIt(g,n)==INVALID) ++n;
83 83
        start=n;
84 84
      }
85 85
      if (start!=INVALID) {
86 86
        for (NodeIt n(g); n!=INVALID; ++n) narc[n]=OutArcIt(g,n);
87 87
        while (narc[start]!=INVALID) {
88 88
          euler.push_back(narc[start]);
89 89
          Node next=g.target(narc[start]);
90 90
          ++narc[start];
91 91
          start=next;
92 92
        }
93 93
      }
94 94
    }
95 95

	
96 96
    ///Arc conversion
97 97
    operator Arc() { return euler.empty()?INVALID:euler.front(); }
98 98
    ///Compare with \c INVALID
99 99
    bool operator==(Invalid) { return euler.empty(); }
100 100
    ///Compare with \c INVALID
101 101
    bool operator!=(Invalid) { return !euler.empty(); }
102 102

	
103 103
    ///Next arc of the tour
104 104

	
105 105
    ///Next arc of the tour
106 106
    ///
107 107
    DiEulerIt &operator++() {
108 108
      Node s=g.target(euler.front());
109 109
      euler.pop_front();
110 110
      typename std::list<Arc>::iterator next=euler.begin();
111 111
      while(narc[s]!=INVALID) {
112 112
        euler.insert(next,narc[s]);
113 113
        Node n=g.target(narc[s]);
114 114
        ++narc[s];
115 115
        s=n;
116 116
      }
117 117
      return *this;
118 118
    }
119 119
    ///Postfix incrementation
120 120

	
121 121
    /// Postfix incrementation.
122 122
    ///
123 123
    ///\warning This incrementation
124 124
    ///returns an \c Arc, not a \ref DiEulerIt, as one may
125 125
    ///expect.
126 126
    Arc operator++(int)
127 127
    {
128 128
      Arc e=*this;
129 129
      ++(*this);
130 130
      return e;
131 131
    }
132 132
  };
133 133

	
134 134
  ///Euler tour iterator for graphs.
135 135

	
136 136
  /// \ingroup graph_properties
137 137
  ///This iterator provides an Euler tour (Eulerian circuit) of an
138 138
  ///\e undirected graph (if there exists) and it converts to the \c Arc
139 139
  ///and \c Edge types of the graph.
140 140
  ///
141 141
  ///For example, if the given graph has an Euler tour (i.e it has only one
142 142
  ///non-trivial component and the degree of each node is even),
143 143
  ///the following code will print the arc IDs according to an
144 144
  ///Euler tour of \c g.
145 145
  ///\code
146 146
  ///  for(EulerIt<ListGraph> e(g); e!=INVALID; ++e) {
147 147
  ///    std::cout << g.id(Edge(e)) << std::eol;
148 148
  ///  }
149 149
  ///\endcode
150 150
  ///Although this iterator is for undirected graphs, it still returns
151 151
  ///arcs in order to indicate the direction of the tour.
152 152
  ///(But arcs convert to edges, of course.)
153 153
  ///
154 154
  ///If \c g has no Euler tour, then the resulted walk will not be closed
155 155
  ///or not contain all edges.
156 156
  template<typename GR>
157 157
  class EulerIt
158 158
  {
159 159
    typedef typename GR::Node Node;
160 160
    typedef typename GR::NodeIt NodeIt;
161 161
    typedef typename GR::Arc Arc;
162 162
    typedef typename GR::Edge Edge;
163 163
    typedef typename GR::ArcIt ArcIt;
164 164
    typedef typename GR::OutArcIt OutArcIt;
165 165
    typedef typename GR::InArcIt InArcIt;
166 166

	
167 167
    const GR &g;
168 168
    typename GR::template NodeMap<OutArcIt> narc;
169 169
    typename GR::template EdgeMap<bool> visited;
170 170
    std::list<Arc> euler;
171 171

	
172 172
  public:
173 173

	
174 174
    ///Constructor
175 175

	
176 176
    ///Constructor.
177 177
    ///\param gr A graph.
178 178
    ///\param start The starting point of the tour. If it is not given,
179 179
    ///the tour will start from the first node that has an incident edge.
180 180
    EulerIt(const GR &gr, typename GR::Node start = INVALID)
181 181
      : g(gr), narc(g), visited(g, false)
182 182
    {
183 183
      if (start==INVALID) {
184 184
        NodeIt n(g);
185 185
        while (n!=INVALID && OutArcIt(g,n)==INVALID) ++n;
186 186
        start=n;
187 187
      }
188 188
      if (start!=INVALID) {
189 189
        for (NodeIt n(g); n!=INVALID; ++n) narc[n]=OutArcIt(g,n);
190 190
        while(narc[start]!=INVALID) {
191 191
          euler.push_back(narc[start]);
192 192
          visited[narc[start]]=true;
193 193
          Node next=g.target(narc[start]);
194 194
          ++narc[start];
195 195
          start=next;
196 196
          while(narc[start]!=INVALID && visited[narc[start]]) ++narc[start];
197 197
        }
198 198
      }
199 199
    }
200 200

	
201 201
    ///Arc conversion
202 202
    operator Arc() const { return euler.empty()?INVALID:euler.front(); }
203 203
    ///Edge conversion
204 204
    operator Edge() const { return euler.empty()?INVALID:euler.front(); }
205 205
    ///Compare with \c INVALID
206 206
    bool operator==(Invalid) const { return euler.empty(); }
207 207
    ///Compare with \c INVALID
208 208
    bool operator!=(Invalid) const { return !euler.empty(); }
209 209

	
210 210
    ///Next arc of the tour
211 211

	
212 212
    ///Next arc of the tour
213 213
    ///
214 214
    EulerIt &operator++() {
215 215
      Node s=g.target(euler.front());
216 216
      euler.pop_front();
217 217
      typename std::list<Arc>::iterator next=euler.begin();
218 218
      while(narc[s]!=INVALID) {
219 219
        while(narc[s]!=INVALID && visited[narc[s]]) ++narc[s];
220 220
        if(narc[s]==INVALID) break;
221 221
        else {
222 222
          euler.insert(next,narc[s]);
223 223
          visited[narc[s]]=true;
224 224
          Node n=g.target(narc[s]);
225 225
          ++narc[s];
226 226
          s=n;
227 227
        }
228 228
      }
229 229
      return *this;
230 230
    }
231 231

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

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

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

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

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

	
34 34
namespace lemon {
35 35

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

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

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

	
82 82
  public:
83 83

	
84 84
    /// \brief Problem type constants for the \c run() function.
85 85
    ///
86 86
    /// Enum type containing the problem type constants that can be
87 87
    /// returned by the \ref run() function of the algorithm.
88 88
    enum ProblemType {
89 89
      /// The problem has no feasible solution (flow).
90 90
      INFEASIBLE,
91 91
      /// The problem has optimal solution (i.e. it is feasible and
92 92
      /// bounded), and the algorithm has found optimal flow and node
93 93
      /// potentials (primal and dual solutions).
94 94
      OPTIMAL,
95 95
      /// The objective function of the problem is unbounded, i.e.
96 96
      /// there is a directed cycle having negative total cost and
97 97
      /// infinite upper bound.
98 98
      UNBOUNDED
99 99
    };
100 100

	
101 101
    /// \brief Constants for selecting the type of the supply constraints.
102 102
    ///
103 103
    /// Enum type containing constants for selecting the supply type,
104 104
    /// i.e. the direction of the inequalities in the supply/demand
105 105
    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
106 106
    ///
107 107
    /// The default supply type is \c GEQ, the \c LEQ type can be
108 108
    /// selected using \ref supplyType().
109 109
    /// The equality form is a special case of both supply types.
110 110
    enum SupplyType {
111 111
      /// This option means that there are <em>"greater or equal"</em>
112 112
      /// supply/demand constraints in the definition of the problem.
113 113
      GEQ,
114 114
      /// This option means that there are <em>"less or equal"</em>
115 115
      /// supply/demand constraints in the definition of the problem.
116 116
      LEQ
117 117
    };
118 118

	
119 119
    /// \brief Constants for selecting the pivot rule.
120 120
    ///
121 121
    /// Enum type containing constants for selecting the pivot rule for
122 122
    /// the \ref run() function.
123 123
    ///
124 124
    /// \ref NetworkSimplex provides five different pivot rule
125 125
    /// implementations that significantly affect the running time
126 126
    /// of the algorithm.
127 127
    /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
128
    /// proved to be the most efficient and the most robust on various
128
    /// turend out to be the most efficient and the most robust on various
129 129
    /// test inputs.
130 130
    /// However, another pivot rule can be selected using the \ref run()
131 131
    /// function with the proper parameter.
132 132
    enum PivotRule {
133 133

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

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

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

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

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

	
162 162
  private:
163 163

	
164 164
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
165 165

	
166 166
    typedef std::vector<int> IntVector;
167 167
    typedef std::vector<Value> ValueVector;
168 168
    typedef std::vector<Cost> CostVector;
169 169
    typedef std::vector<signed char> CharVector;
170
    // Note: vector<signed char> is used instead of vector<ArcState> and 
170
    // Note: vector<signed char> is used instead of vector<ArcState> and
171 171
    // vector<ArcDirection> for efficiency reasons
172 172

	
173 173
    // State constants for arcs
174 174
    enum ArcState {
175 175
      STATE_UPPER = -1,
176 176
      STATE_TREE  =  0,
177 177
      STATE_LOWER =  1
178 178
    };
179 179

	
180 180
    // Direction constants for tree arcs
181 181
    enum ArcDirection {
182 182
      DIR_DOWN = -1,
183 183
      DIR_UP   =  1
184 184
    };
185 185

	
186 186
  private:
187 187

	
188 188
    // Data related to the underlying digraph
189 189
    const GR &_graph;
190 190
    int _node_num;
191 191
    int _arc_num;
192 192
    int _all_arc_num;
193 193
    int _search_arc_num;
194 194

	
195 195
    // Parameters of the problem
196 196
    bool _have_lower;
197 197
    SupplyType _stype;
198 198
    Value _sum_supply;
199 199

	
200 200
    // Data structures for storing the digraph
201 201
    IntNodeMap _node_id;
202 202
    IntArcMap _arc_id;
203 203
    IntVector _source;
204 204
    IntVector _target;
205 205
    bool _arc_mixing;
206 206

	
207 207
    // Node and arc data
208 208
    ValueVector _lower;
209 209
    ValueVector _upper;
210 210
    ValueVector _cap;
211 211
    CostVector _cost;
212 212
    ValueVector _supply;
213 213
    ValueVector _flow;
214 214
    CostVector _pi;
215 215

	
216 216
    // Data for storing the spanning tree structure
217 217
    IntVector _parent;
218 218
    IntVector _pred;
219 219
    IntVector _thread;
220 220
    IntVector _rev_thread;
221 221
    IntVector _succ_num;
222 222
    IntVector _last_succ;
223 223
    CharVector _pred_dir;
224 224
    CharVector _state;
225 225
    IntVector _dirty_revs;
226 226
    int _root;
227 227

	
228 228
    // Temporary data used in the current pivot iteration
229 229
    int in_arc, join, u_in, v_in, u_out, v_out;
230 230
    Value delta;
231 231

	
232 232
    const Value MAX;
233 233

	
234 234
  public:
235 235

	
236 236
    /// \brief Constant for infinite upper bounds (capacities).
237 237
    ///
238 238
    /// Constant for infinite upper bounds (capacities).
239 239
    /// It is \c std::numeric_limits<Value>::infinity() if available,
240 240
    /// \c std::numeric_limits<Value>::max() otherwise.
241 241
    const Value INF;
242 242

	
243 243
  private:
244 244

	
245 245
    // Implementation of the First Eligible pivot rule
246 246
    class FirstEligiblePivotRule
247 247
    {
248 248
    private:
249 249

	
250 250
      // References to the NetworkSimplex class
251 251
      const IntVector  &_source;
252 252
      const IntVector  &_target;
253 253
      const CostVector &_cost;
254 254
      const CharVector &_state;
255 255
      const CostVector &_pi;
256 256
      int &_in_arc;
257 257
      int _search_arc_num;
258 258

	
259 259
      // Pivot rule data
260 260
      int _next_arc;
261 261

	
262 262
    public:
263 263

	
264 264
      // Constructor
265 265
      FirstEligiblePivotRule(NetworkSimplex &ns) :
266 266
        _source(ns._source), _target(ns._target),
267 267
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
268 268
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
269 269
        _next_arc(0)
270 270
      {}
271 271

	
272 272
      // Find next entering arc
273 273
      bool findEnteringArc() {
274 274
        Cost c;
275 275
        for (int e = _next_arc; e != _search_arc_num; ++e) {
276 276
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
277 277
          if (c < 0) {
278 278
            _in_arc = e;
279 279
            _next_arc = e + 1;
280 280
            return true;
281 281
          }
282 282
        }
283 283
        for (int e = 0; e != _next_arc; ++e) {
284 284
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
285 285
          if (c < 0) {
286 286
            _in_arc = e;
287 287
            _next_arc = e + 1;
288 288
            return true;
289 289
          }
290 290
        }
291 291
        return false;
292 292
      }
293 293

	
294 294
    }; //class FirstEligiblePivotRule
295 295

	
296 296

	
297 297
    // Implementation of the Best Eligible pivot rule
298 298
    class BestEligiblePivotRule
299 299
    {
300 300
    private:
301 301

	
302 302
      // References to the NetworkSimplex class
303 303
      const IntVector  &_source;
304 304
      const IntVector  &_target;
305 305
      const CostVector &_cost;
306 306
      const CharVector &_state;
307 307
      const CostVector &_pi;
308 308
      int &_in_arc;
309 309
      int _search_arc_num;
310 310

	
311 311
    public:
312 312

	
313 313
      // Constructor
314 314
      BestEligiblePivotRule(NetworkSimplex &ns) :
315 315
        _source(ns._source), _target(ns._target),
316 316
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
317 317
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
318 318
      {}
319 319

	
320 320
      // Find next entering arc
321 321
      bool findEnteringArc() {
322 322
        Cost c, min = 0;
323 323
        for (int e = 0; e != _search_arc_num; ++e) {
324 324
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
325 325
          if (c < min) {
326 326
            min = c;
327 327
            _in_arc = e;
328 328
          }
329 329
        }
330 330
        return min < 0;
331 331
      }
332 332

	
333 333
    }; //class BestEligiblePivotRule
334 334

	
335 335

	
336 336
    // Implementation of the Block Search pivot rule
337 337
    class BlockSearchPivotRule
338 338
    {
339 339
    private:
340 340

	
341 341
      // References to the NetworkSimplex class
342 342
      const IntVector  &_source;
343 343
      const IntVector  &_target;
344 344
      const CostVector &_cost;
345 345
      const CharVector &_state;
346 346
      const CostVector &_pi;
347 347
      int &_in_arc;
348 348
      int _search_arc_num;
349 349

	
350 350
      // Pivot rule data
351 351
      int _block_size;
352 352
      int _next_arc;
353 353

	
354 354
    public:
355 355

	
356 356
      // Constructor
357 357
      BlockSearchPivotRule(NetworkSimplex &ns) :
358 358
        _source(ns._source), _target(ns._target),
359 359
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
360 360
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
361 361
        _next_arc(0)
362 362
      {
... ...
@@ -545,401 +545,403 @@
545 545

	
546 546
    public:
547 547

	
548 548
      // Constructor
549 549
      AlteringListPivotRule(NetworkSimplex &ns) :
550 550
        _source(ns._source), _target(ns._target),
551 551
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
552 552
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
553 553
        _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
554 554
      {
555 555
        // The main parameters of the pivot rule
556 556
        const double BLOCK_SIZE_FACTOR = 1.0;
557 557
        const int MIN_BLOCK_SIZE = 10;
558 558
        const double HEAD_LENGTH_FACTOR = 0.1;
559 559
        const int MIN_HEAD_LENGTH = 3;
560 560

	
561 561
        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
562 562
                                    std::sqrt(double(_search_arc_num))),
563 563
                                MIN_BLOCK_SIZE );
564 564
        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
565 565
                                 MIN_HEAD_LENGTH );
566 566
        _candidates.resize(_head_length + _block_size);
567 567
        _curr_length = 0;
568 568
      }
569 569

	
570 570
      // Find next entering arc
571 571
      bool findEnteringArc() {
572 572
        // Check the current candidate list
573 573
        int e;
574 574
        Cost c;
575 575
        for (int i = 0; i != _curr_length; ++i) {
576 576
          e = _candidates[i];
577 577
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
578 578
          if (c < 0) {
579 579
            _cand_cost[e] = c;
580 580
          } else {
581 581
            _candidates[i--] = _candidates[--_curr_length];
582 582
          }
583 583
        }
584 584

	
585 585
        // Extend the list
586 586
        int cnt = _block_size;
587 587
        int limit = _head_length;
588 588

	
589 589
        for (e = _next_arc; e != _search_arc_num; ++e) {
590 590
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
591 591
          if (c < 0) {
592 592
            _cand_cost[e] = c;
593 593
            _candidates[_curr_length++] = e;
594 594
          }
595 595
          if (--cnt == 0) {
596 596
            if (_curr_length > limit) goto search_end;
597 597
            limit = 0;
598 598
            cnt = _block_size;
599 599
          }
600 600
        }
601 601
        for (e = 0; e != _next_arc; ++e) {
602 602
          _cand_cost[e] = _state[e] *
603 603
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
604 604
          if (_cand_cost[e] < 0) {
605 605
            _candidates[_curr_length++] = e;
606 606
          }
607 607
          if (--cnt == 0) {
608 608
            if (_curr_length > limit) goto search_end;
609 609
            limit = 0;
610 610
            cnt = _block_size;
611 611
          }
612 612
        }
613 613
        if (_curr_length == 0) return false;
614 614

	
615 615
      search_end:
616 616

	
617 617
        // Make heap of the candidate list (approximating a partial sort)
618 618
        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
619 619
                   _sort_func );
620 620

	
621 621
        // Pop the first element of the heap
622 622
        _in_arc = _candidates[0];
623 623
        _next_arc = e;
624 624
        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
625 625
                  _sort_func );
626 626
        _curr_length = std::min(_head_length, _curr_length - 1);
627 627
        return true;
628 628
      }
629 629

	
630 630
    }; //class AlteringListPivotRule
631 631

	
632 632
  public:
633 633

	
634 634
    /// \brief Constructor.
635 635
    ///
636 636
    /// The constructor of the class.
637 637
    ///
638 638
    /// \param graph The digraph the algorithm runs on.
639 639
    /// \param arc_mixing Indicate if the arcs will be stored in a
640 640
    /// mixed order in the internal data structure.
641 641
    /// In general, it leads to similar performance as using the original
642 642
    /// arc order, but it makes the algorithm more robust and in special
643 643
    /// cases, even significantly faster. Therefore, it is enabled by default.
644 644
    NetworkSimplex(const GR& graph, bool arc_mixing = true) :
645 645
      _graph(graph), _node_id(graph), _arc_id(graph),
646 646
      _arc_mixing(arc_mixing),
647 647
      MAX(std::numeric_limits<Value>::max()),
648 648
      INF(std::numeric_limits<Value>::has_infinity ?
649 649
          std::numeric_limits<Value>::infinity() : MAX)
650 650
    {
651 651
      // Check the number types
652 652
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
653 653
        "The flow type of NetworkSimplex must be signed");
654 654
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
655 655
        "The cost type of NetworkSimplex must be signed");
656 656

	
657 657
      // Reset data structures
658 658
      reset();
659 659
    }
660 660

	
661 661
    /// \name Parameters
662 662
    /// The parameters of the algorithm can be specified using these
663 663
    /// functions.
664 664

	
665 665
    /// @{
666 666

	
667 667
    /// \brief Set the lower bounds on the arcs.
668 668
    ///
669 669
    /// This function sets the lower bounds on the arcs.
670 670
    /// If it is not used before calling \ref run(), the lower bounds
671 671
    /// will be set to zero on all arcs.
672 672
    ///
673 673
    /// \param map An arc map storing the lower bounds.
674 674
    /// Its \c Value type must be convertible to the \c Value type
675 675
    /// of the algorithm.
676 676
    ///
677 677
    /// \return <tt>(*this)</tt>
678 678
    template <typename LowerMap>
679 679
    NetworkSimplex& lowerMap(const LowerMap& map) {
680 680
      _have_lower = true;
681 681
      for (ArcIt a(_graph); a != INVALID; ++a) {
682 682
        _lower[_arc_id[a]] = map[a];
683 683
      }
684 684
      return *this;
685 685
    }
686 686

	
687 687
    /// \brief Set the upper bounds (capacities) on the arcs.
688 688
    ///
689 689
    /// This function sets the upper bounds (capacities) on the arcs.
690 690
    /// If it is not used before calling \ref run(), the upper bounds
691 691
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
692 692
    /// unbounded from above).
693 693
    ///
694 694
    /// \param map An arc map storing the upper bounds.
695 695
    /// Its \c Value type must be convertible to the \c Value type
696 696
    /// of the algorithm.
697 697
    ///
698 698
    /// \return <tt>(*this)</tt>
699 699
    template<typename UpperMap>
700 700
    NetworkSimplex& upperMap(const UpperMap& map) {
701 701
      for (ArcIt a(_graph); a != INVALID; ++a) {
702 702
        _upper[_arc_id[a]] = map[a];
703 703
      }
704 704
      return *this;
705 705
    }
706 706

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

	
726 726
    /// \brief Set the supply values of the nodes.
727 727
    ///
728 728
    /// This function sets the supply values of the nodes.
729 729
    /// If neither this function nor \ref stSupply() is used before
730 730
    /// calling \ref run(), the supply of each node will be set to zero.
731 731
    ///
732 732
    /// \param map A node map storing the supply values.
733 733
    /// Its \c Value type must be convertible to the \c Value type
734 734
    /// of the algorithm.
735 735
    ///
736 736
    /// \return <tt>(*this)</tt>
737
    ///
738
    /// \sa supplyType()
737 739
    template<typename SupplyMap>
738 740
    NetworkSimplex& supplyMap(const SupplyMap& map) {
739 741
      for (NodeIt n(_graph); n != INVALID; ++n) {
740 742
        _supply[_node_id[n]] = map[n];
741 743
      }
742 744
      return *this;
743 745
    }
744 746

	
745 747
    /// \brief Set single source and target nodes and a supply value.
746 748
    ///
747 749
    /// This function sets a single source node and a single target node
748 750
    /// and the required flow value.
749 751
    /// If neither this function nor \ref supplyMap() is used before
750 752
    /// calling \ref run(), the supply of each node will be set to zero.
751 753
    ///
752 754
    /// Using this function has the same effect as using \ref supplyMap()
753
    /// with such a map in which \c k is assigned to \c s, \c -k is
755
    /// with a map in which \c k is assigned to \c s, \c -k is
754 756
    /// assigned to \c t and all other nodes have zero supply value.
755 757
    ///
756 758
    /// \param s The source node.
757 759
    /// \param t The target node.
758 760
    /// \param k The required amount of flow from node \c s to node \c t
759 761
    /// (i.e. the supply of \c s and the demand of \c t).
760 762
    ///
761 763
    /// \return <tt>(*this)</tt>
762 764
    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
763 765
      for (int i = 0; i != _node_num; ++i) {
764 766
        _supply[i] = 0;
765 767
      }
766 768
      _supply[_node_id[s]] =  k;
767 769
      _supply[_node_id[t]] = -k;
768 770
      return *this;
769 771
    }
770 772

	
771 773
    /// \brief Set the type of the supply constraints.
772 774
    ///
773 775
    /// This function sets the type of the supply/demand constraints.
774 776
    /// If it is not used before calling \ref run(), the \ref GEQ supply
775 777
    /// type will be used.
776 778
    ///
777 779
    /// For more information, see \ref SupplyType.
778 780
    ///
779 781
    /// \return <tt>(*this)</tt>
780 782
    NetworkSimplex& supplyType(SupplyType supply_type) {
781 783
      _stype = supply_type;
782 784
      return *this;
783 785
    }
784 786

	
785 787
    /// @}
786 788

	
787 789
    /// \name Execution Control
788 790
    /// The algorithm can be executed using \ref run().
789 791

	
790 792
    /// @{
791 793

	
792 794
    /// \brief Run the algorithm.
793 795
    ///
794 796
    /// This function runs the algorithm.
795 797
    /// The paramters can be specified using functions \ref lowerMap(),
796 798
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
797 799
    /// \ref supplyType().
798 800
    /// For example,
799 801
    /// \code
800 802
    ///   NetworkSimplex<ListDigraph> ns(graph);
801 803
    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
802 804
    ///     .supplyMap(sup).run();
803 805
    /// \endcode
804 806
    ///
805 807
    /// This function can be called more than once. All the given parameters
806 808
    /// are kept for the next call, unless \ref resetParams() or \ref reset()
807 809
    /// is used, thus only the modified parameters have to be set again.
808 810
    /// If the underlying digraph was also modified after the construction
809 811
    /// of the class (or the last \ref reset() call), then the \ref reset()
810 812
    /// function must be called.
811 813
    ///
812 814
    /// \param pivot_rule The pivot rule that will be used during the
813 815
    /// algorithm. For more information, see \ref PivotRule.
814 816
    ///
815 817
    /// \return \c INFEASIBLE if no feasible flow exists,
816 818
    /// \n \c OPTIMAL if the problem has optimal solution
817 819
    /// (i.e. it is feasible and bounded), and the algorithm has found
818 820
    /// optimal flow and node potentials (primal and dual solutions),
819 821
    /// \n \c UNBOUNDED if the objective function of the problem is
820 822
    /// unbounded, i.e. there is a directed cycle having negative total
821 823
    /// cost and infinite upper bound.
822 824
    ///
823 825
    /// \see ProblemType, PivotRule
824 826
    /// \see resetParams(), reset()
825 827
    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
826 828
      if (!init()) return INFEASIBLE;
827 829
      return start(pivot_rule);
828 830
    }
829 831

	
830 832
    /// \brief Reset all the parameters that have been given before.
831 833
    ///
832 834
    /// This function resets all the paramaters that have been given
833 835
    /// before using functions \ref lowerMap(), \ref upperMap(),
834 836
    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
835 837
    ///
836 838
    /// It is useful for multiple \ref run() calls. Basically, all the given
837 839
    /// parameters are kept for the next \ref run() call, unless
838 840
    /// \ref resetParams() or \ref reset() is used.
839 841
    /// If the underlying digraph was also modified after the construction
840 842
    /// of the class or the last \ref reset() call, then the \ref reset()
841 843
    /// function must be used, otherwise \ref resetParams() is sufficient.
842 844
    ///
843 845
    /// For example,
844 846
    /// \code
845 847
    ///   NetworkSimplex<ListDigraph> ns(graph);
846 848
    ///
847 849
    ///   // First run
848 850
    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
849 851
    ///     .supplyMap(sup).run();
850 852
    ///
851 853
    ///   // Run again with modified cost map (resetParams() is not called,
852 854
    ///   // so only the cost map have to be set again)
853 855
    ///   cost[e] += 100;
854 856
    ///   ns.costMap(cost).run();
855 857
    ///
856 858
    ///   // Run again from scratch using resetParams()
857 859
    ///   // (the lower bounds will be set to zero on all arcs)
858 860
    ///   ns.resetParams();
859 861
    ///   ns.upperMap(capacity).costMap(cost)
860 862
    ///     .supplyMap(sup).run();
861 863
    /// \endcode
862 864
    ///
863 865
    /// \return <tt>(*this)</tt>
864 866
    ///
865 867
    /// \see reset(), run()
866 868
    NetworkSimplex& resetParams() {
867 869
      for (int i = 0; i != _node_num; ++i) {
868 870
        _supply[i] = 0;
869 871
      }
870 872
      for (int i = 0; i != _arc_num; ++i) {
871 873
        _lower[i] = 0;
872 874
        _upper[i] = INF;
873 875
        _cost[i] = 1;
874 876
      }
875 877
      _have_lower = false;
876 878
      _stype = GEQ;
877 879
      return *this;
878 880
    }
879 881

	
880 882
    /// \brief Reset the internal data structures and all the parameters
881 883
    /// that have been given before.
882 884
    ///
883 885
    /// This function resets the internal data structures and all the
884 886
    /// paramaters that have been given before using functions \ref lowerMap(),
885 887
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
886 888
    /// \ref supplyType().
887 889
    ///
888 890
    /// It is useful for multiple \ref run() calls. Basically, all the given
889 891
    /// parameters are kept for the next \ref run() call, unless
890 892
    /// \ref resetParams() or \ref reset() is used.
891 893
    /// If the underlying digraph was also modified after the construction
892 894
    /// of the class or the last \ref reset() call, then the \ref reset()
893 895
    /// function must be used, otherwise \ref resetParams() is sufficient.
894 896
    ///
895 897
    /// See \ref resetParams() for examples.
896 898
    ///
897 899
    /// \return <tt>(*this)</tt>
898 900
    ///
899 901
    /// \see resetParams(), run()
900 902
    NetworkSimplex& reset() {
901 903
      // Resize vectors
902 904
      _node_num = countNodes(_graph);
903 905
      _arc_num = countArcs(_graph);
904 906
      int all_node_num = _node_num + 1;
905 907
      int max_arc_num = _arc_num + 2 * _node_num;
906 908

	
907 909
      _source.resize(max_arc_num);
908 910
      _target.resize(max_arc_num);
909 911

	
910 912
      _lower.resize(_arc_num);
911 913
      _upper.resize(_arc_num);
912 914
      _cap.resize(max_arc_num);
913 915
      _cost.resize(max_arc_num);
914 916
      _supply.resize(all_node_num);
915 917
      _flow.resize(max_arc_num);
916 918
      _pi.resize(all_node_num);
917 919

	
918 920
      _parent.resize(all_node_num);
919 921
      _pred.resize(all_node_num);
920 922
      _pred_dir.resize(all_node_num);
921 923
      _thread.resize(all_node_num);
922 924
      _rev_thread.resize(all_node_num);
923 925
      _succ_num.resize(all_node_num);
924 926
      _last_succ.resize(all_node_num);
925 927
      _state.resize(max_arc_num);
926 928

	
927 929
      // Copy the graph
928 930
      int i = 0;
929 931
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
930 932
        _node_id[n] = i;
931 933
      }
932 934
      if (_arc_mixing) {
933 935
        // Store the arcs in a mixed order
934 936
        const int skip = std::max(_arc_num / _node_num, 3);
935 937
        int i = 0, j = 0;
936 938
        for (ArcIt a(_graph); a != INVALID; ++a) {
937 939
          _arc_id[a] = i;
938 940
          _source[i] = _node_id[_graph.source(a)];
939 941
          _target[i] = _node_id[_graph.target(a)];
940 942
          if ((i += skip) >= _arc_num) i = ++j;
941 943
        }
942 944
      } else {
943 945
        // Store the arcs in the original order
944 946
        int i = 0;
945 947
        for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
0 comments (0 inline)