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
... ...
@@ -25,742 +25,742 @@
25 25

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	
265 265
/**
266 266
@defgroup 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.
735 735
*/
736 736

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

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

	
745 745
/**
746 746
@defgroup tools Standalone Utility Applications
747 747

	
748 748
Some utility applications are listed here.
749 749

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

	
754 754
/**
755 755
\anchor demoprograms
756 756

	
757 757
@defgroup demos Demo Programs
758 758

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

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

	
766 766
}
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];
618 618
        int bi = _arc_idb[a];
619 619
        _reverse[fi] = bi;
620 620
        _reverse[bi] = fi;
621 621
      }
622 622

	
623 623
      // Reset parameters
624 624
      resetParams();
625 625
      return *this;
626 626
    }
627 627

	
628 628
    /// @}
629 629

	
630 630
    /// \name Query Functions
631 631
    /// The results of the algorithm can be obtained using these
632 632
    /// functions.\n
633 633
    /// The \ref run() function must be called before using them.
634 634

	
635 635
    /// @{
636 636

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

	
663 663
#ifndef DOXYGEN
664 664
    Cost totalCost() const {
665 665
      return totalCost<Cost>();
666 666
    }
667 667
#endif
668 668

	
669 669
    /// \brief Return the flow on the given arc.
670 670
    ///
671 671
    /// This function returns the flow on the given arc.
672 672
    ///
673 673
    /// \pre \ref run() must be called before using this function.
674 674
    Value flow(const Arc& a) const {
675 675
      return _res_cap[_arc_idb[a]];
676 676
    }
677 677

	
678 678
    /// \brief Return the flow map (the primal solution).
679 679
    ///
680 680
    /// This function copies the flow value on each arc into the given
681 681
    /// map. The \c Value type of the algorithm must be convertible to
682 682
    /// the \c Value type of the map.
683 683
    ///
684 684
    /// \pre \ref run() must be called before using this function.
685 685
    template <typename FlowMap>
686 686
    void flowMap(FlowMap &map) const {
687 687
      for (ArcIt a(_graph); a != INVALID; ++a) {
688 688
        map.set(a, _res_cap[_arc_idb[a]]);
689 689
      }
690 690
    }
691 691

	
692 692
    /// \brief Return the potential (dual value) of the given node.
693 693
    ///
694 694
    /// This function returns the potential (dual value) of the
695 695
    /// given node.
696 696
    ///
697 697
    /// \pre \ref run() must be called before using this function.
698 698
    Cost potential(const Node& n) const {
699 699
      return _pi[_node_id[n]];
700 700
    }
701 701

	
702 702
    /// \brief Return the potential map (the dual solution).
703 703
    ///
704 704
    /// This function copies the potential (dual value) of each node
705 705
    /// into the given map.
706 706
    /// The \c Cost type of the algorithm must be convertible to the
707 707
    /// \c Value type of the map.
708 708
    ///
709 709
    /// \pre \ref run() must be called before using this function.
710 710
    template <typename PotentialMap>
711 711
    void potentialMap(PotentialMap &map) const {
712 712
      for (NodeIt n(_graph); n != INVALID; ++n) {
713 713
        map.set(n, _pi[_node_id[n]]);
714 714
      }
715 715
    }
716 716

	
717 717
    /// @}
718 718

	
719 719
  private:
720 720

	
721 721
    // Initialize the algorithm
722 722
    ProblemType init() {
723 723
      if (_node_num <= 1) return INFEASIBLE;
724 724

	
725 725
      // Check the sum of supply values
726 726
      _sum_supply = 0;
727 727
      for (int i = 0; i != _root; ++i) {
728 728
        _sum_supply += _supply[i];
729 729
      }
730 730
      if (_sum_supply > 0) return INFEASIBLE;
731 731

	
732 732
      // Initialize vectors
733 733
      for (int i = 0; i != _root; ++i) {
734 734
        _pi[i] = 0;
735 735
        _excess[i] = _supply[i];
736 736
      }
737 737

	
738 738
      // Remove non-zero lower bounds
739 739
      const Value MAX = std::numeric_limits<Value>::max();
740 740
      int last_out;
741 741
      if (_have_lower) {
742 742
        for (int i = 0; i != _root; ++i) {
743 743
          last_out = _first_out[i+1];
744 744
          for (int j = _first_out[i]; j != last_out; ++j) {
745 745
            if (_forward[j]) {
746 746
              Value c = _lower[j];
747 747
              if (c >= 0) {
748 748
                _res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF;
749 749
              } else {
750 750
                _res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF;
751 751
              }
752 752
              _excess[i] -= c;
753 753
              _excess[_target[j]] += c;
754 754
            } else {
755 755
              _res_cap[j] = 0;
756 756
            }
757 757
          }
758 758
        }
759 759
      } else {
760 760
        for (int j = 0; j != _res_arc_num; ++j) {
761 761
          _res_cap[j] = _forward[j] ? _upper[j] : 0;
762 762
        }
763 763
      }
764 764

	
765 765
      // Handle negative costs
766 766
      for (int i = 0; i != _root; ++i) {
767 767
        last_out = _first_out[i+1] - 1;
768 768
        for (int j = _first_out[i]; j != last_out; ++j) {
769 769
          Value rc = _res_cap[j];
770 770
          if (_cost[j] < 0 && rc > 0) {
771 771
            if (rc >= MAX) return UNBOUNDED;
772 772
            _excess[i] -= rc;
773 773
            _excess[_target[j]] += rc;
774 774
            _res_cap[j] = 0;
775 775
            _res_cap[_reverse[j]] += rc;
776 776
          }
777 777
        }
778 778
      }
779 779

	
780 780
      // Handle GEQ supply type
781 781
      if (_sum_supply < 0) {
782 782
        _pi[_root] = 0;
783 783
        _excess[_root] = -_sum_supply;
784 784
        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
785 785
          int ra = _reverse[a];
786 786
          _res_cap[a] = -_sum_supply + 1;
787 787
          _res_cap[ra] = 0;
788 788
          _cost[a] = 0;
789 789
          _cost[ra] = 0;
790 790
        }
791 791
      } else {
792 792
        _pi[_root] = 0;
793 793
        _excess[_root] = 0;
794 794
        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
795 795
          int ra = _reverse[a];
796 796
          _res_cap[a] = 1;
797 797
          _res_cap[ra] = 0;
798 798
          _cost[a] = 0;
799 799
          _cost[ra] = 0;
800 800
        }
801 801
      }
802 802

	
803 803
      // Initialize delta value
804 804
      if (_factor > 1) {
805 805
        // With scaling
806 806
        Value max_sup = 0, max_dem = 0, max_cap = 0;
807 807
        for (int i = 0; i != _root; ++i) {
808 808
          Value ex = _excess[i];
809 809
          if ( ex > max_sup) max_sup =  ex;
Ignore white space 6 line context
... ...
@@ -66,769 +66,769 @@
66 66
#else
67 67
  extern const Invalid INVALID;
68 68
#endif
69 69

	
70 70
  /// \addtogroup gutils
71 71
  /// @{
72 72

	
73 73
  ///Create convenience typedefs for the digraph types and iterators
74 74

	
75 75
  ///This \c \#define creates convenient type definitions for the following
76 76
  ///types of \c Digraph: \c Node,  \c NodeIt, \c Arc, \c ArcIt, \c InArcIt,
77 77
  ///\c OutArcIt, \c BoolNodeMap, \c IntNodeMap, \c DoubleNodeMap,
78 78
  ///\c BoolArcMap, \c IntArcMap, \c DoubleArcMap.
79 79
  ///
80 80
  ///\note If the graph type is a dependent type, ie. the graph type depend
81 81
  ///on a template parameter, then use \c TEMPLATE_DIGRAPH_TYPEDEFS()
82 82
  ///macro.
83 83
#define DIGRAPH_TYPEDEFS(Digraph)                                       \
84 84
  typedef Digraph::Node Node;                                           \
85 85
  typedef Digraph::NodeIt NodeIt;                                       \
86 86
  typedef Digraph::Arc Arc;                                             \
87 87
  typedef Digraph::ArcIt ArcIt;                                         \
88 88
  typedef Digraph::InArcIt InArcIt;                                     \
89 89
  typedef Digraph::OutArcIt OutArcIt;                                   \
90 90
  typedef Digraph::NodeMap<bool> BoolNodeMap;                           \
91 91
  typedef Digraph::NodeMap<int> IntNodeMap;                             \
92 92
  typedef Digraph::NodeMap<double> DoubleNodeMap;                       \
93 93
  typedef Digraph::ArcMap<bool> BoolArcMap;                             \
94 94
  typedef Digraph::ArcMap<int> IntArcMap;                               \
95 95
  typedef Digraph::ArcMap<double> DoubleArcMap
96 96

	
97 97
  ///Create convenience typedefs for the digraph types and iterators
98 98

	
99 99
  ///\see DIGRAPH_TYPEDEFS
100 100
  ///
101 101
  ///\note Use this macro, if the graph type is a dependent type,
102 102
  ///ie. the graph type depend on a template parameter.
103 103
#define TEMPLATE_DIGRAPH_TYPEDEFS(Digraph)                              \
104 104
  typedef typename Digraph::Node Node;                                  \
105 105
  typedef typename Digraph::NodeIt NodeIt;                              \
106 106
  typedef typename Digraph::Arc Arc;                                    \
107 107
  typedef typename Digraph::ArcIt ArcIt;                                \
108 108
  typedef typename Digraph::InArcIt InArcIt;                            \
109 109
  typedef typename Digraph::OutArcIt OutArcIt;                          \
110 110
  typedef typename Digraph::template NodeMap<bool> BoolNodeMap;         \
111 111
  typedef typename Digraph::template NodeMap<int> IntNodeMap;           \
112 112
  typedef typename Digraph::template NodeMap<double> DoubleNodeMap;     \
113 113
  typedef typename Digraph::template ArcMap<bool> BoolArcMap;           \
114 114
  typedef typename Digraph::template ArcMap<int> IntArcMap;             \
115 115
  typedef typename Digraph::template ArcMap<double> DoubleArcMap
116 116

	
117 117
  ///Create convenience typedefs for the graph types and iterators
118 118

	
119 119
  ///This \c \#define creates the same convenient type definitions as defined
120 120
  ///by \ref DIGRAPH_TYPEDEFS(Graph) and six more, namely it creates
121 121
  ///\c Edge, \c EdgeIt, \c IncEdgeIt, \c BoolEdgeMap, \c IntEdgeMap,
122 122
  ///\c DoubleEdgeMap.
123 123
  ///
124 124
  ///\note If the graph type is a dependent type, ie. the graph type depend
125 125
  ///on a template parameter, then use \c TEMPLATE_GRAPH_TYPEDEFS()
126 126
  ///macro.
127 127
#define GRAPH_TYPEDEFS(Graph)                                           \
128 128
  DIGRAPH_TYPEDEFS(Graph);                                              \
129 129
  typedef Graph::Edge Edge;                                             \
130 130
  typedef Graph::EdgeIt EdgeIt;                                         \
131 131
  typedef Graph::IncEdgeIt IncEdgeIt;                                   \
132 132
  typedef Graph::EdgeMap<bool> BoolEdgeMap;                             \
133 133
  typedef Graph::EdgeMap<int> IntEdgeMap;                               \
134 134
  typedef Graph::EdgeMap<double> DoubleEdgeMap
135 135

	
136 136
  ///Create convenience typedefs for the graph types and iterators
137 137

	
138 138
  ///\see GRAPH_TYPEDEFS
139 139
  ///
140 140
  ///\note Use this macro, if the graph type is a dependent type,
141 141
  ///ie. the graph type depend on a template parameter.
142 142
#define TEMPLATE_GRAPH_TYPEDEFS(Graph)                                  \
143 143
  TEMPLATE_DIGRAPH_TYPEDEFS(Graph);                                     \
144 144
  typedef typename Graph::Edge Edge;                                    \
145 145
  typedef typename Graph::EdgeIt EdgeIt;                                \
146 146
  typedef typename Graph::IncEdgeIt IncEdgeIt;                          \
147 147
  typedef typename Graph::template EdgeMap<bool> BoolEdgeMap;           \
148 148
  typedef typename Graph::template EdgeMap<int> IntEdgeMap;             \
149 149
  typedef typename Graph::template EdgeMap<double> DoubleEdgeMap
150 150

	
151 151
  /// \brief Function to count the items in a graph.
152 152
  ///
153 153
  /// This function counts the items (nodes, arcs etc.) in a graph.
154 154
  /// The complexity of the function is linear because
155 155
  /// it iterates on all of the items.
156 156
  template <typename Graph, typename Item>
157 157
  inline int countItems(const Graph& g) {
158 158
    typedef typename ItemSetTraits<Graph, Item>::ItemIt ItemIt;
159 159
    int num = 0;
160 160
    for (ItemIt it(g); it != INVALID; ++it) {
161 161
      ++num;
162 162
    }
163 163
    return num;
164 164
  }
165 165

	
166 166
  // Node counting:
167 167

	
168 168
  namespace _core_bits {
169 169

	
170 170
    template <typename Graph, typename Enable = void>
171 171
    struct CountNodesSelector {
172 172
      static int count(const Graph &g) {
173 173
        return countItems<Graph, typename Graph::Node>(g);
174 174
      }
175 175
    };
176 176

	
177 177
    template <typename Graph>
178 178
    struct CountNodesSelector<
179 179
      Graph, typename
180 180
      enable_if<typename Graph::NodeNumTag, void>::type>
181 181
    {
182 182
      static int count(const Graph &g) {
183 183
        return g.nodeNum();
184 184
      }
185 185
    };
186 186
  }
187 187

	
188 188
  /// \brief Function to count the nodes in the graph.
189 189
  ///
190 190
  /// This function counts the nodes in the graph.
191 191
  /// The complexity of the function is <em>O</em>(<em>n</em>), but for some
192 192
  /// graph structures it is specialized to run in <em>O</em>(1).
193 193
  ///
194 194
  /// \note If the graph contains a \c nodeNum() member function and a
195 195
  /// \c NodeNumTag tag then this function calls directly the member
196 196
  /// function to query the cardinality of the node set.
197 197
  template <typename Graph>
198 198
  inline int countNodes(const Graph& g) {
199 199
    return _core_bits::CountNodesSelector<Graph>::count(g);
200 200
  }
201 201

	
202 202
  // Arc counting:
203 203

	
204 204
  namespace _core_bits {
205 205

	
206 206
    template <typename Graph, typename Enable = void>
207 207
    struct CountArcsSelector {
208 208
      static int count(const Graph &g) {
209 209
        return countItems<Graph, typename Graph::Arc>(g);
210 210
      }
211 211
    };
212 212

	
213 213
    template <typename Graph>
214 214
    struct CountArcsSelector<
215 215
      Graph,
216 216
      typename enable_if<typename Graph::ArcNumTag, void>::type>
217 217
    {
218 218
      static int count(const Graph &g) {
219 219
        return g.arcNum();
220 220
      }
221 221
    };
222 222
  }
223 223

	
224 224
  /// \brief Function to count the arcs in the graph.
225 225
  ///
226 226
  /// This function counts the arcs in the graph.
227 227
  /// The complexity of the function is <em>O</em>(<em>m</em>), but for some
228 228
  /// graph structures it is specialized to run in <em>O</em>(1).
229 229
  ///
230 230
  /// \note If the graph contains a \c arcNum() member function and a
231 231
  /// \c ArcNumTag tag then this function calls directly the member
232 232
  /// function to query the cardinality of the arc set.
233 233
  template <typename Graph>
234 234
  inline int countArcs(const Graph& g) {
235 235
    return _core_bits::CountArcsSelector<Graph>::count(g);
236 236
  }
237 237

	
238 238
  // Edge counting:
239 239

	
240 240
  namespace _core_bits {
241 241

	
242 242
    template <typename Graph, typename Enable = void>
243 243
    struct CountEdgesSelector {
244 244
      static int count(const Graph &g) {
245 245
        return countItems<Graph, typename Graph::Edge>(g);
246 246
      }
247 247
    };
248 248

	
249 249
    template <typename Graph>
250 250
    struct CountEdgesSelector<
251 251
      Graph,
252 252
      typename enable_if<typename Graph::EdgeNumTag, void>::type>
253 253
    {
254 254
      static int count(const Graph &g) {
255 255
        return g.edgeNum();
256 256
      }
257 257
    };
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.
643 643
    void run() {
644 644
      NodeRefMap nodeRefMap(_from);
645 645
      ArcRefMap arcRefMap(_from);
646 646
      _core_bits::DigraphCopySelector<To>::
647 647
        copy(_from, _to, nodeRefMap, arcRefMap);
648 648
      for (int i = 0; i < int(_node_maps.size()); ++i) {
649 649
        _node_maps[i]->copy(_from, nodeRefMap);
650 650
      }
651 651
      for (int i = 0; i < int(_arc_maps.size()); ++i) {
652 652
        _arc_maps[i]->copy(_from, arcRefMap);
653 653
      }
654 654
    }
655 655

	
656 656
  protected:
657 657

	
658 658
    const From& _from;
659 659
    To& _to;
660 660

	
661 661
    std::vector<_core_bits::MapCopyBase<From, Node, NodeRefMap>* >
662 662
      _node_maps;
663 663

	
664 664
    std::vector<_core_bits::MapCopyBase<From, Arc, ArcRefMap>* >
665 665
      _arc_maps;
666 666

	
667 667
  };
668 668

	
669 669
  /// \brief Copy a digraph to another digraph.
670 670
  ///
671 671
  /// This function copies a digraph to another digraph.
672 672
  /// The complete usage of it is detailed in the DigraphCopy class, but
673 673
  /// a short example shows a basic work:
674 674
  ///\code
675 675
  /// digraphCopy(src, trg).nodeRef(nr).arcCrossRef(acr).run();
676 676
  ///\endcode
677 677
  ///
678 678
  /// After the copy the \c nr map will contain the mapping from the
679 679
  /// nodes of the \c from digraph to the nodes of the \c to digraph and
680 680
  /// \c acr will contain the mapping from the arcs of the \c to digraph
681 681
  /// to the arcs of the \c from digraph.
682 682
  ///
683 683
  /// \see DigraphCopy
684 684
  template <typename From, typename To>
685 685
  DigraphCopy<From, To> digraphCopy(const From& from, To& to) {
686 686
    return DigraphCopy<From, To>(from, to);
687 687
  }
688 688

	
689 689
  /// \brief Class to copy a graph.
690 690
  ///
691 691
  /// Class to copy a graph to another graph (duplicate a graph). The
692 692
  /// simplest way of using it is through the \c graphCopy() function.
693 693
  ///
694 694
  /// This class not only make a copy of a graph, but it can create
695 695
  /// references and cross references between the nodes, edges and arcs of
696 696
  /// the two graphs, and it can copy maps for using with the newly created
697 697
  /// graph.
698 698
  ///
699 699
  /// To make a copy from a graph, first an instance of GraphCopy
700 700
  /// should be created, then the data belongs to the graph should
701 701
  /// assigned to copy. In the end, the \c run() member should be
702 702
  /// called.
703 703
  ///
704 704
  /// The next code copies a graph with several data:
705 705
  ///\code
706 706
  ///  GraphCopy<OrigGraph, NewGraph> cg(orig_graph, new_graph);
707 707
  ///  // Create references for the nodes
708 708
  ///  OrigGraph::NodeMap<NewGraph::Node> nr(orig_graph);
709 709
  ///  cg.nodeRef(nr);
710 710
  ///  // Create cross references (inverse) for the edges
711 711
  ///  NewGraph::EdgeMap<OrigGraph::Edge> ecr(new_graph);
712 712
  ///  cg.edgeCrossRef(ecr);
713 713
  ///  // Copy an edge map
714 714
  ///  OrigGraph::EdgeMap<double> oemap(orig_graph);
715 715
  ///  NewGraph::EdgeMap<double> nemap(new_graph);
716 716
  ///  cg.edgeMap(oemap, nemap);
717 717
  ///  // Copy a node
718 718
  ///  OrigGraph::Node on;
719 719
  ///  NewGraph::Node nn;
720 720
  ///  cg.node(on, nn);
721 721
  ///  // Execute copying
722 722
  ///  cg.run();
723 723
  ///\endcode
724 724
  template <typename From, typename To>
725 725
  class GraphCopy {
726 726
  private:
727 727

	
728 728
    typedef typename From::Node Node;
729 729
    typedef typename From::NodeIt NodeIt;
730 730
    typedef typename From::Arc Arc;
731 731
    typedef typename From::ArcIt ArcIt;
732 732
    typedef typename From::Edge Edge;
733 733
    typedef typename From::EdgeIt EdgeIt;
734 734

	
735 735
    typedef typename To::Node TNode;
736 736
    typedef typename To::Arc TArc;
737 737
    typedef typename To::Edge TEdge;
738 738

	
739 739
    typedef typename From::template NodeMap<TNode> NodeRefMap;
740 740
    typedef typename From::template EdgeMap<TEdge> EdgeRefMap;
741 741

	
742 742
    struct ArcRefMap {
743 743
      ArcRefMap(const From& from, const To& to,
744 744
                const EdgeRefMap& edge_ref, const NodeRefMap& node_ref)
745 745
        : _from(from), _to(to),
746 746
          _edge_ref(edge_ref), _node_ref(node_ref) {}
747 747

	
748 748
      typedef typename From::Arc Key;
749 749
      typedef typename To::Arc Value;
750 750

	
751 751
      Value operator[](const Key& key) const {
752 752
        bool forward = _from.u(key) != _from.v(key) ?
753 753
          _node_ref[_from.source(key)] ==
754 754
          _to.source(_to.direct(_edge_ref[key], true)) :
755 755
          _from.direction(key);
756 756
        return _to.direct(_edge_ref[key], forward);
757 757
      }
758 758

	
759 759
      const From& _from;
760 760
      const To& _to;
761 761
      const EdgeRefMap& _edge_ref;
762 762
      const NodeRefMap& _node_ref;
763 763
    };
764 764

	
765 765
  public:
766 766

	
767 767
    /// \brief Constructor of GraphCopy.
768 768
    ///
769 769
    /// Constructor of GraphCopy for copying the content of the
770 770
    /// \c from graph into the \c to graph.
771 771
    GraphCopy(const From& from, To& to)
772 772
      : _from(from), _to(to) {}
773 773

	
774 774
    /// \brief Destructor of GraphCopy
775 775
    ///
776 776
    /// Destructor of GraphCopy.
777 777
    ~GraphCopy() {
778 778
      for (int i = 0; i < int(_node_maps.size()); ++i) {
779 779
        delete _node_maps[i];
780 780
      }
781 781
      for (int i = 0; i < int(_arc_maps.size()); ++i) {
782 782
        delete _arc_maps[i];
783 783
      }
784 784
      for (int i = 0; i < int(_edge_maps.size()); ++i) {
785 785
        delete _edge_maps[i];
786 786
      }
787 787
    }
788 788

	
789 789
    /// \brief Copy the node references into the given map.
790 790
    ///
791 791
    /// This function copies the node references into the given map.
792 792
    /// The parameter should be a map, whose key type is the Node type of
793 793
    /// the source graph, while the value type is the Node type of the
794 794
    /// destination graph.
795 795
    template <typename NodeRef>
796 796
    GraphCopy& nodeRef(NodeRef& map) {
797 797
      _node_maps.push_back(new _core_bits::RefCopy<From, Node,
798 798
                           NodeRefMap, NodeRef>(map));
799 799
      return *this;
800 800
    }
801 801

	
802 802
    /// \brief Copy the node cross references into the given map.
803 803
    ///
804 804
    /// This function copies the node cross references (reverse references)
805 805
    /// into the given map. The parameter should be a map, whose key type
806 806
    /// is the Node type of the destination graph, while the value type is
807 807
    /// the Node type of the source graph.
808 808
    template <typename NodeCrossRef>
809 809
    GraphCopy& nodeCrossRef(NodeCrossRef& map) {
810 810
      _node_maps.push_back(new _core_bits::CrossRefCopy<From, Node,
811 811
                           NodeRefMap, NodeCrossRef>(map));
812 812
      return *this;
813 813
    }
814 814

	
815 815
    /// \brief Make a copy of the given node map.
816 816
    ///
817 817
    /// This function makes a copy of the given node map for the newly
818 818
    /// created graph.
819 819
    /// The key type of the new map \c tmap should be the Node type of the
820 820
    /// destination graph, and the key type of the original map \c map
821 821
    /// should be the Node type of the source graph.
822 822
    template <typename FromMap, typename ToMap>
823 823
    GraphCopy& nodeMap(const FromMap& map, ToMap& tmap) {
824 824
      _node_maps.push_back(new _core_bits::MapCopy<From, Node,
825 825
                           NodeRefMap, FromMap, ToMap>(map, tmap));
826 826
      return *this;
827 827
    }
828 828

	
829 829
    /// \brief Make a copy of the given node.
830 830
    ///
831 831
    /// This function makes a copy of the given node.
832 832
    GraphCopy& node(const Node& node, TNode& tnode) {
833 833
      _node_maps.push_back(new _core_bits::ItemCopy<From, Node,
834 834
                           NodeRefMap, TNode>(node, tnode));
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_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;
643 646
      }
644 647
      _first_out[i] = j;
645 648
      _first_out[_res_node_num] = k;
646 649
      for (ArcIt a(_graph); a != INVALID; ++a) {
647 650
        int fi = _arc_idf[a];
648 651
        int bi = _arc_idb[a];
649 652
        _reverse[fi] = bi;
650 653
        _reverse[bi] = fi;
651 654
      }
652 655

	
653 656
      // Reset parameters
654 657
      resetParams();
655 658
      return *this;
656 659
    }
657 660

	
658 661
    /// @}
659 662

	
660 663
    /// \name Query Functions
661 664
    /// The results of the algorithm can be obtained using these
662 665
    /// functions.\n
663 666
    /// The \ref run() function must be called before using them.
664 667

	
665 668
    /// @{
666 669

	
667 670
    /// \brief Return the total cost of the found flow.
668 671
    ///
669 672
    /// This function returns the total cost of the found flow.
670 673
    /// Its complexity is O(e).
671 674
    ///
672 675
    /// \note The return type of the function can be specified as a
673 676
    /// template parameter. For example,
674 677
    /// \code
675 678
    ///   cs.totalCost<double>();
676 679
    /// \endcode
677 680
    /// It is useful if the total cost cannot be stored in the \c Cost
678 681
    /// type of the algorithm, which is the default return type of the
679 682
    /// function.
680 683
    ///
681 684
    /// \pre \ref run() must be called before using this function.
682 685
    template <typename Number>
683 686
    Number totalCost() const {
684 687
      Number c = 0;
685 688
      for (ArcIt a(_graph); a != INVALID; ++a) {
686 689
        int i = _arc_idb[a];
687 690
        c += static_cast<Number>(_res_cap[i]) *
688 691
             (-static_cast<Number>(_scost[i]));
689 692
      }
690 693
      return c;
691 694
    }
692 695

	
693 696
#ifndef DOXYGEN
694 697
    Cost totalCost() const {
695 698
      return totalCost<Cost>();
696 699
    }
697 700
#endif
698 701

	
699 702
    /// \brief Return the flow on the given arc.
700 703
    ///
701 704
    /// This function returns the flow on the given arc.
702 705
    ///
703 706
    /// \pre \ref run() must be called before using this function.
704 707
    Value flow(const Arc& a) const {
705 708
      return _res_cap[_arc_idb[a]];
706 709
    }
707 710

	
708 711
    /// \brief Return the flow map (the primal solution).
709 712
    ///
710 713
    /// This function copies the flow value on each arc into the given
711 714
    /// map. The \c Value type of the algorithm must be convertible to
712 715
    /// the \c Value type of the map.
713 716
    ///
714 717
    /// \pre \ref run() must be called before using this function.
715 718
    template <typename FlowMap>
716 719
    void flowMap(FlowMap &map) const {
717 720
      for (ArcIt a(_graph); a != INVALID; ++a) {
718 721
        map.set(a, _res_cap[_arc_idb[a]]);
719 722
      }
720 723
    }
721 724

	
722 725
    /// \brief Return the potential (dual value) of the given node.
723 726
    ///
724 727
    /// This function returns the potential (dual value) of the
725 728
    /// given node.
726 729
    ///
727 730
    /// \pre \ref run() must be called before using this function.
728 731
    Cost potential(const Node& n) const {
729 732
      return static_cast<Cost>(_pi[_node_id[n]]);
730 733
    }
731 734

	
732 735
    /// \brief Return the potential map (the dual solution).
733 736
    ///
734 737
    /// This function copies the potential (dual value) of each node
735 738
    /// into the given map.
736 739
    /// The \c Cost type of the algorithm must be convertible to the
737 740
    /// \c Value type of the map.
738 741
    ///
739 742
    /// \pre \ref run() must be called before using this function.
740 743
    template <typename PotentialMap>
741 744
    void potentialMap(PotentialMap &map) const {
742 745
      for (NodeIt n(_graph); n != INVALID; ++n) {
743 746
        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
744 747
      }
745 748
    }
746 749

	
747 750
    /// @}
748 751

	
749 752
  private:
750 753

	
751 754
    // Initialize the algorithm
752 755
    ProblemType init() {
753 756
      if (_res_node_num <= 1) return INFEASIBLE;
754 757

	
755 758
      // Check the sum of supply values
756 759
      _sum_supply = 0;
757 760
      for (int i = 0; i != _root; ++i) {
758 761
        _sum_supply += _supply[i];
759 762
      }
760 763
      if (_sum_supply > 0) return INFEASIBLE;
761 764

	
762 765

	
763 766
      // Initialize vectors
764 767
      for (int i = 0; i != _res_node_num; ++i) {
765 768
        _pi[i] = 0;
766 769
        _excess[i] = _supply[i];
767 770
      }
768 771

	
769 772
      // Remove infinite upper bounds and check negative arcs
770 773
      const Value MAX = std::numeric_limits<Value>::max();
771 774
      int last_out;
772 775
      if (_have_lower) {
773 776
        for (int i = 0; i != _root; ++i) {
774 777
          last_out = _first_out[i+1];
775 778
          for (int j = _first_out[i]; j != last_out; ++j) {
776 779
            if (_forward[j]) {
777 780
              Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
778 781
              if (c >= MAX) return UNBOUNDED;
779 782
              _excess[i] -= c;
780 783
              _excess[_target[j]] += c;
781 784
            }
782 785
          }
783 786
        }
784 787
      } else {
785 788
        for (int i = 0; i != _root; ++i) {
786 789
          last_out = _first_out[i+1];
787 790
          for (int j = _first_out[i]; j != last_out; ++j) {
788 791
            if (_forward[j] && _scost[j] < 0) {
789 792
              Value c = _upper[j];
790 793
              if (c >= MAX) return UNBOUNDED;
791 794
              _excess[i] -= c;
792 795
              _excess[_target[j]] += c;
793 796
            }
794 797
          }
795 798
        }
796 799
      }
797 800
      Value ex, max_cap = 0;
798 801
      for (int i = 0; i != _res_node_num; ++i) {
799 802
        ex = _excess[i];
800 803
        _excess[i] = 0;
801 804
        if (ex < 0) max_cap -= ex;
802 805
      }
803 806
      for (int j = 0; j != _res_arc_num; ++j) {
804 807
        if (_upper[j] >= MAX) _upper[j] = max_cap;
805 808
      }
806 809

	
807 810
      // Initialize the large cost vector and the epsilon parameter
808 811
      _epsilon = 0;
809 812
      LargeCost lc;
810 813
      for (int i = 0; i != _root; ++i) {
811 814
        last_out = _first_out[i+1];
812 815
        for (int j = _first_out[i]; j != last_out; ++j) {
813 816
          lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
814 817
          _cost[j] = lc;
815 818
          if (lc > _epsilon) _epsilon = lc;
816 819
        }
817 820
      }
818 821
      _epsilon /= _alpha;
819 822

	
820 823
      // Initialize maps for Circulation and remove non-zero lower bounds
821 824
      ConstMap<Arc, Value> low(0);
822 825
      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
823 826
      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
824 827
      ValueArcMap cap(_graph), flow(_graph);
825 828
      ValueNodeMap sup(_graph);
826 829
      for (NodeIt n(_graph); n != INVALID; ++n) {
827 830
        sup[n] = _supply[_node_id[n]];
828 831
      }
829 832
      if (_have_lower) {
830 833
        for (ArcIt a(_graph); a != INVALID; ++a) {
831 834
          int j = _arc_idf[a];
832 835
          Value c = _lower[j];
833 836
          cap[a] = _upper[j] - c;
834 837
          sup[_graph.source(a)] -= c;
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;
545 544
        _target[k] = i;
546 545
        _reverse[k] = j;
547 546
        ++j; ++k;
548 547
      }
549 548
      _first_out[i] = j;
550 549
      _first_out[_res_node_num] = k;
551 550
      for (ArcIt a(_graph); a != INVALID; ++a) {
552 551
        int fi = _arc_idf[a];
553 552
        int bi = _arc_idb[a];
554 553
        _reverse[fi] = bi;
555 554
        _reverse[bi] = fi;
556 555
      }
557 556

	
558 557
      // Reset parameters
559 558
      resetParams();
560 559
      return *this;
561 560
    }
562 561

	
563 562
    /// @}
564 563

	
565 564
    /// \name Query Functions
566 565
    /// The results of the algorithm can be obtained using these
567 566
    /// functions.\n
568 567
    /// The \ref run() function must be called before using them.
569 568

	
570 569
    /// @{
571 570

	
572 571
    /// \brief Return the total cost of the found flow.
573 572
    ///
574 573
    /// This function returns the total cost of the found flow.
575 574
    /// Its complexity is O(e).
576 575
    ///
577 576
    /// \note The return type of the function can be specified as a
578 577
    /// template parameter. For example,
579 578
    /// \code
580 579
    ///   cc.totalCost<double>();
581 580
    /// \endcode
582 581
    /// It is useful if the total cost cannot be stored in the \c Cost
583 582
    /// type of the algorithm, which is the default return type of the
584 583
    /// function.
585 584
    ///
586 585
    /// \pre \ref run() must be called before using this function.
587 586
    template <typename Number>
588 587
    Number totalCost() const {
589 588
      Number c = 0;
590 589
      for (ArcIt a(_graph); a != INVALID; ++a) {
591 590
        int i = _arc_idb[a];
592 591
        c += static_cast<Number>(_res_cap[i]) *
593 592
             (-static_cast<Number>(_cost[i]));
594 593
      }
595 594
      return c;
596 595
    }
597 596

	
598 597
#ifndef DOXYGEN
599 598
    Cost totalCost() const {
600 599
      return totalCost<Cost>();
601 600
    }
602 601
#endif
603 602

	
604 603
    /// \brief Return the flow on the given arc.
605 604
    ///
606 605
    /// This function returns the flow on the given arc.
607 606
    ///
608 607
    /// \pre \ref run() must be called before using this function.
609 608
    Value flow(const Arc& a) const {
610 609
      return _res_cap[_arc_idb[a]];
611 610
    }
612 611

	
613 612
    /// \brief Return the flow map (the primal solution).
614 613
    ///
615 614
    /// This function copies the flow value on each arc into the given
616 615
    /// map. The \c Value type of the algorithm must be convertible to
617 616
    /// the \c Value type of the map.
618 617
    ///
619 618
    /// \pre \ref run() must be called before using this function.
620 619
    template <typename FlowMap>
621 620
    void flowMap(FlowMap &map) const {
622 621
      for (ArcIt a(_graph); a != INVALID; ++a) {
623 622
        map.set(a, _res_cap[_arc_idb[a]]);
624 623
      }
625 624
    }
626 625

	
627 626
    /// \brief Return the potential (dual value) of the given node.
628 627
    ///
629 628
    /// This function returns the potential (dual value) of the
630 629
    /// given node.
631 630
    ///
632 631
    /// \pre \ref run() must be called before using this function.
633 632
    Cost potential(const Node& n) const {
634 633
      return static_cast<Cost>(_pi[_node_id[n]]);
635 634
    }
636 635

	
637 636
    /// \brief Return the potential map (the dual solution).
638 637
    ///
639 638
    /// This function copies the potential (dual value) of each node
640 639
    /// into the given map.
641 640
    /// The \c Cost type of the algorithm must be convertible to the
642 641
    /// \c Value type of the map.
643 642
    ///
644 643
    /// \pre \ref run() must be called before using this function.
645 644
    template <typename PotentialMap>
646 645
    void potentialMap(PotentialMap &map) const {
647 646
      for (NodeIt n(_graph); n != INVALID; ++n) {
648 647
        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
649 648
      }
650 649
    }
651 650

	
652 651
    /// @}
653 652

	
654 653
  private:
655 654

	
656 655
    // Initialize the algorithm
657 656
    ProblemType init() {
658 657
      if (_res_node_num <= 1) return INFEASIBLE;
659 658

	
660 659
      // Check the sum of supply values
661 660
      _sum_supply = 0;
662 661
      for (int i = 0; i != _root; ++i) {
663 662
        _sum_supply += _supply[i];
664 663
      }
665 664
      if (_sum_supply > 0) return INFEASIBLE;
666 665

	
667 666

	
668 667
      // Initialize vectors
669 668
      for (int i = 0; i != _res_node_num; ++i) {
670 669
        _pi[i] = 0;
671 670
      }
672 671
      ValueVector excess(_supply);
673 672

	
674 673
      // Remove infinite upper bounds and check negative arcs
675 674
      const Value MAX = std::numeric_limits<Value>::max();
676 675
      int last_out;
677 676
      if (_have_lower) {
678 677
        for (int i = 0; i != _root; ++i) {
679 678
          last_out = _first_out[i+1];
680 679
          for (int j = _first_out[i]; j != last_out; ++j) {
681 680
            if (_forward[j]) {
682 681
              Value c = _cost[j] < 0 ? _upper[j] : _lower[j];
683 682
              if (c >= MAX) return UNBOUNDED;
684 683
              excess[i] -= c;
685 684
              excess[_target[j]] += c;
686 685
            }
687 686
          }
688 687
        }
689 688
      } else {
690 689
        for (int i = 0; i != _root; ++i) {
691 690
          last_out = _first_out[i+1];
692 691
          for (int j = _first_out[i]; j != last_out; ++j) {
693 692
            if (_forward[j] && _cost[j] < 0) {
694 693
              Value c = _upper[j];
695 694
              if (c >= MAX) return UNBOUNDED;
696 695
              excess[i] -= c;
697 696
              excess[_target[j]] += c;
698 697
            }
699 698
          }
700 699
        }
701 700
      }
702 701
      Value ex, max_cap = 0;
703 702
      for (int i = 0; i != _res_node_num; ++i) {
704 703
        ex = excess[i];
705 704
        if (ex < 0) max_cap -= ex;
706 705
      }
707 706
      for (int j = 0; j != _res_arc_num; ++j) {
708 707
        if (_upper[j] >= MAX) _upper[j] = max_cap;
709 708
      }
710 709

	
711 710
      // Initialize maps for Circulation and remove non-zero lower bounds
712 711
      ConstMap<Arc, Value> low(0);
713 712
      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
714 713
      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
715 714
      ValueArcMap cap(_graph), flow(_graph);
716 715
      ValueNodeMap sup(_graph);
717 716
      for (NodeIt n(_graph); n != INVALID; ++n) {
718 717
        sup[n] = _supply[_node_id[n]];
719 718
      }
720 719
      if (_have_lower) {
721 720
        for (ArcIt a(_graph); a != INVALID; ++a) {
722 721
          int j = _arc_idf[a];
723 722
          Value c = _lower[j];
724 723
          cap[a] = _upper[j] - c;
725 724
          sup[_graph.source(a)] -= c;
726 725
          sup[_graph.target(a)] += c;
727 726
        }
728 727
      } else {
729 728
        for (ArcIt a(_graph); a != INVALID; ++a) {
730 729
          cap[a] = _upper[_arc_idf[a]];
731 730
        }
732 731
      }
733 732

	
734 733
      // Find a feasible flow using Circulation
735 734
      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
736 735
        circ(_graph, low, cap, sup);
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

	
232 232
    ///Postfix incrementation
233 233

	
234 234
    /// Postfix incrementation.
235 235
    ///
236 236
    ///\warning This incrementation returns an \c Arc (which converts to
237 237
    ///an \c Edge), not an \ref EulerIt, as one may expect.
238 238
    Arc operator++(int)
239 239
    {
240 240
      Arc e=*this;
241 241
      ++(*this);
242 242
      return e;
243 243
    }
244 244
  };
245 245

	
246 246

	
247 247
  ///Check if the given graph is Eulerian
248 248

	
249 249
  /// \ingroup graph_properties
250 250
  ///This function checks if the given graph is Eulerian.
251 251
  ///It works for both directed and undirected graphs.
252 252
  ///
253 253
  ///By definition, a digraph is called \e Eulerian if
254 254
  ///and only if it is connected and the number of incoming and outgoing
255 255
  ///arcs are the same for each node.
256 256
  ///Similarly, an undirected graph is called \e Eulerian if
257 257
  ///and only if it is connected and the number of incident edges is even
258 258
  ///for each node.
259 259
  ///
260 260
  ///\note There are (di)graphs that are not Eulerian, but still have an
261 261
  /// Euler tour, since they may contain isolated nodes.
262 262
  ///
263 263
  ///\sa DiEulerIt, EulerIt
264 264
  template<typename GR>
265 265
#ifdef DOXYGEN
266 266
  bool
267 267
#else
268 268
  typename enable_if<UndirectedTagIndicator<GR>,bool>::type
269 269
  eulerian(const GR &g)
270 270
  {
271 271
    for(typename GR::NodeIt n(g);n!=INVALID;++n)
272 272
      if(countIncEdges(g,n)%2) return false;
273 273
    return connected(g);
274 274
  }
275 275
  template<class GR>
276 276
  typename disable_if<UndirectedTagIndicator<GR>,bool>::type
277 277
#endif
278 278
  eulerian(const GR &g)
279 279
  {
280 280
    for(typename GR::NodeIt n(g);n!=INVALID;++n)
281 281
      if(countInArcs(g,n)!=countOutArcs(g,n)) return false;
282 282
    return connected(undirector(g));
283 283
  }
284 284

	
285 285
}
286 286

	
287 287
#endif
Ignore white space 768 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
      {
363 363
        // The main parameters of the pivot rule
364 364
        const double BLOCK_SIZE_FACTOR = 1.0;
365 365
        const int MIN_BLOCK_SIZE = 10;
366 366

	
367 367
        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
368 368
                                    std::sqrt(double(_search_arc_num))),
369 369
                                MIN_BLOCK_SIZE );
370 370
      }
371 371

	
372 372
      // Find next entering arc
373 373
      bool findEnteringArc() {
374 374
        Cost c, min = 0;
375 375
        int cnt = _block_size;
376 376
        int e;
377 377
        for (e = _next_arc; e != _search_arc_num; ++e) {
378 378
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
379 379
          if (c < min) {
380 380
            min = c;
381 381
            _in_arc = e;
382 382
          }
383 383
          if (--cnt == 0) {
384 384
            if (min < 0) goto search_end;
385 385
            cnt = _block_size;
386 386
          }
387 387
        }
388 388
        for (e = 0; e != _next_arc; ++e) {
389 389
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
390 390
          if (c < min) {
391 391
            min = c;
392 392
            _in_arc = e;
393 393
          }
394 394
          if (--cnt == 0) {
395 395
            if (min < 0) goto search_end;
396 396
            cnt = _block_size;
397 397
          }
398 398
        }
399 399
        if (min >= 0) return false;
400 400

	
401 401
      search_end:
402 402
        _next_arc = e;
403 403
        return true;
404 404
      }
405 405

	
406 406
    }; //class BlockSearchPivotRule
407 407

	
408 408

	
409 409
    // Implementation of the Candidate List pivot rule
410 410
    class CandidateListPivotRule
411 411
    {
412 412
    private:
413 413

	
414 414
      // References to the NetworkSimplex class
415 415
      const IntVector  &_source;
416 416
      const IntVector  &_target;
417 417
      const CostVector &_cost;
418 418
      const CharVector &_state;
419 419
      const CostVector &_pi;
420 420
      int &_in_arc;
421 421
      int _search_arc_num;
422 422

	
423 423
      // Pivot rule data
424 424
      IntVector _candidates;
425 425
      int _list_length, _minor_limit;
426 426
      int _curr_length, _minor_count;
427 427
      int _next_arc;
428 428

	
429 429
    public:
430 430

	
431 431
      /// Constructor
432 432
      CandidateListPivotRule(NetworkSimplex &ns) :
433 433
        _source(ns._source), _target(ns._target),
434 434
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
435 435
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
436 436
        _next_arc(0)
437 437
      {
438 438
        // The main parameters of the pivot rule
439 439
        const double LIST_LENGTH_FACTOR = 0.25;
440 440
        const int MIN_LIST_LENGTH = 10;
441 441
        const double MINOR_LIMIT_FACTOR = 0.1;
442 442
        const int MIN_MINOR_LIMIT = 3;
443 443

	
444 444
        _list_length = std::max( int(LIST_LENGTH_FACTOR *
445 445
                                     std::sqrt(double(_search_arc_num))),
446 446
                                 MIN_LIST_LENGTH );
447 447
        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
448 448
                                 MIN_MINOR_LIMIT );
449 449
        _curr_length = _minor_count = 0;
450 450
        _candidates.resize(_list_length);
451 451
      }
452 452

	
453 453
      /// Find next entering arc
454 454
      bool findEnteringArc() {
455 455
        Cost min, c;
456 456
        int e;
457 457
        if (_curr_length > 0 && _minor_count < _minor_limit) {
458 458
          // Minor iteration: select the best eligible arc from the
459 459
          // current candidate list
460 460
          ++_minor_count;
461 461
          min = 0;
462 462
          for (int i = 0; i < _curr_length; ++i) {
463 463
            e = _candidates[i];
464 464
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
465 465
            if (c < min) {
466 466
              min = c;
467 467
              _in_arc = e;
468 468
            }
469 469
            else if (c >= 0) {
470 470
              _candidates[i--] = _candidates[--_curr_length];
471 471
            }
472 472
          }
473 473
          if (min < 0) return true;
474 474
        }
475 475

	
476 476
        // Major iteration: build a new candidate list
477 477
        min = 0;
478 478
        _curr_length = 0;
479 479
        for (e = _next_arc; e != _search_arc_num; ++e) {
480 480
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
481 481
          if (c < 0) {
482 482
            _candidates[_curr_length++] = e;
483 483
            if (c < min) {
484 484
              min = c;
485 485
              _in_arc = e;
486 486
            }
487 487
            if (_curr_length == _list_length) goto search_end;
488 488
          }
489 489
        }
490 490
        for (e = 0; e != _next_arc; ++e) {
491 491
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
492 492
          if (c < 0) {
493 493
            _candidates[_curr_length++] = e;
494 494
            if (c < min) {
495 495
              min = c;
496 496
              _in_arc = e;
497 497
            }
498 498
            if (_curr_length == _list_length) goto search_end;
499 499
          }
500 500
        }
501 501
        if (_curr_length == 0) return false;
502 502

	
503 503
      search_end:
504 504
        _minor_count = 1;
505 505
        _next_arc = e;
506 506
        return true;
507 507
      }
508 508

	
509 509
    }; //class CandidateListPivotRule
510 510

	
511 511

	
512 512
    // Implementation of the Altering Candidate List pivot rule
513 513
    class AlteringListPivotRule
514 514
    {
515 515
    private:
516 516

	
517 517
      // References to the NetworkSimplex class
518 518
      const IntVector  &_source;
519 519
      const IntVector  &_target;
520 520
      const CostVector &_cost;
521 521
      const CharVector &_state;
522 522
      const CostVector &_pi;
523 523
      int &_in_arc;
524 524
      int _search_arc_num;
525 525

	
526 526
      // Pivot rule data
527 527
      int _block_size, _head_length, _curr_length;
528 528
      int _next_arc;
529 529
      IntVector _candidates;
530 530
      CostVector _cand_cost;
531 531

	
532 532
      // Functor class to compare arcs during sort of the candidate list
533 533
      class SortFunc
534 534
      {
535 535
      private:
536 536
        const CostVector &_map;
537 537
      public:
538 538
        SortFunc(const CostVector &map) : _map(map) {}
539 539
        bool operator()(int left, int right) {
540 540
          return _map[left] > _map[right];
541 541
        }
542 542
      };
543 543

	
544 544
      SortFunc _sort_func;
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) {
946 948
          _arc_id[a] = i;
947 949
          _source[i] = _node_id[_graph.source(a)];
948 950
          _target[i] = _node_id[_graph.target(a)];
949 951
        }
950 952
      }
951 953

	
952 954
      // Reset parameters
953 955
      resetParams();
954 956
      return *this;
955 957
    }
956 958

	
957 959
    /// @}
958 960

	
959 961
    /// \name Query Functions
960 962
    /// The results of the algorithm can be obtained using these
961 963
    /// functions.\n
962 964
    /// The \ref run() function must be called before using them.
963 965

	
964 966
    /// @{
965 967

	
966 968
    /// \brief Return the total cost of the found flow.
967 969
    ///
968 970
    /// This function returns the total cost of the found flow.
969 971
    /// Its complexity is O(e).
970 972
    ///
971 973
    /// \note The return type of the function can be specified as a
972 974
    /// template parameter. For example,
973 975
    /// \code
974 976
    ///   ns.totalCost<double>();
975 977
    /// \endcode
976 978
    /// It is useful if the total cost cannot be stored in the \c Cost
977 979
    /// type of the algorithm, which is the default return type of the
978 980
    /// function.
979 981
    ///
980 982
    /// \pre \ref run() must be called before using this function.
981 983
    template <typename Number>
982 984
    Number totalCost() const {
983 985
      Number c = 0;
984 986
      for (ArcIt a(_graph); a != INVALID; ++a) {
985 987
        int i = _arc_id[a];
986 988
        c += Number(_flow[i]) * Number(_cost[i]);
987 989
      }
988 990
      return c;
989 991
    }
990 992

	
991 993
#ifndef DOXYGEN
992 994
    Cost totalCost() const {
993 995
      return totalCost<Cost>();
994 996
    }
995 997
#endif
996 998

	
997 999
    /// \brief Return the flow on the given arc.
998 1000
    ///
999 1001
    /// This function returns the flow on the given arc.
1000 1002
    ///
1001 1003
    /// \pre \ref run() must be called before using this function.
1002 1004
    Value flow(const Arc& a) const {
1003 1005
      return _flow[_arc_id[a]];
1004 1006
    }
1005 1007

	
1006 1008
    /// \brief Return the flow map (the primal solution).
1007 1009
    ///
1008 1010
    /// This function copies the flow value on each arc into the given
1009 1011
    /// map. The \c Value type of the algorithm must be convertible to
1010 1012
    /// the \c Value type of the map.
1011 1013
    ///
1012 1014
    /// \pre \ref run() must be called before using this function.
1013 1015
    template <typename FlowMap>
1014 1016
    void flowMap(FlowMap &map) const {
1015 1017
      for (ArcIt a(_graph); a != INVALID; ++a) {
1016 1018
        map.set(a, _flow[_arc_id[a]]);
1017 1019
      }
1018 1020
    }
1019 1021

	
1020 1022
    /// \brief Return the potential (dual value) of the given node.
1021 1023
    ///
1022 1024
    /// This function returns the potential (dual value) of the
1023 1025
    /// given node.
1024 1026
    ///
1025 1027
    /// \pre \ref run() must be called before using this function.
1026 1028
    Cost potential(const Node& n) const {
1027 1029
      return _pi[_node_id[n]];
1028 1030
    }
1029 1031

	
1030 1032
    /// \brief Return the potential map (the dual solution).
1031 1033
    ///
1032 1034
    /// This function copies the potential (dual value) of each node
1033 1035
    /// into the given map.
1034 1036
    /// The \c Cost type of the algorithm must be convertible to the
1035 1037
    /// \c Value type of the map.
1036 1038
    ///
1037 1039
    /// \pre \ref run() must be called before using this function.
1038 1040
    template <typename PotentialMap>
1039 1041
    void potentialMap(PotentialMap &map) const {
1040 1042
      for (NodeIt n(_graph); n != INVALID; ++n) {
1041 1043
        map.set(n, _pi[_node_id[n]]);
1042 1044
      }
1043 1045
    }
1044 1046

	
1045 1047
    /// @}
1046 1048

	
1047 1049
  private:
1048 1050

	
1049 1051
    // Initialize internal data structures
1050 1052
    bool init() {
1051 1053
      if (_node_num == 0) return false;
1052 1054

	
1053 1055
      // Check the sum of supply values
1054 1056
      _sum_supply = 0;
1055 1057
      for (int i = 0; i != _node_num; ++i) {
1056 1058
        _sum_supply += _supply[i];
1057 1059
      }
1058 1060
      if ( !((_stype == GEQ && _sum_supply <= 0) ||
1059 1061
             (_stype == LEQ && _sum_supply >= 0)) ) return false;
1060 1062

	
1061 1063
      // Remove non-zero lower bounds
1062 1064
      if (_have_lower) {
1063 1065
        for (int i = 0; i != _arc_num; ++i) {
1064 1066
          Value c = _lower[i];
1065 1067
          if (c >= 0) {
1066 1068
            _cap[i] = _upper[i] < MAX ? _upper[i] - c : INF;
1067 1069
          } else {
1068 1070
            _cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF;
1069 1071
          }
1070 1072
          _supply[_source[i]] -= c;
1071 1073
          _supply[_target[i]] += c;
1072 1074
        }
1073 1075
      } else {
1074 1076
        for (int i = 0; i != _arc_num; ++i) {
1075 1077
          _cap[i] = _upper[i];
1076 1078
        }
1077 1079
      }
1078 1080

	
1079 1081
      // Initialize artifical cost
1080 1082
      Cost ART_COST;
1081 1083
      if (std::numeric_limits<Cost>::is_exact) {
1082 1084
        ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
1083 1085
      } else {
1084 1086
        ART_COST = 0;
1085 1087
        for (int i = 0; i != _arc_num; ++i) {
1086 1088
          if (_cost[i] > ART_COST) ART_COST = _cost[i];
1087 1089
        }
1088 1090
        ART_COST = (ART_COST + 1) * _node_num;
1089 1091
      }
1090 1092

	
1091 1093
      // Initialize arc maps
1092 1094
      for (int i = 0; i != _arc_num; ++i) {
1093 1095
        _flow[i] = 0;
1094 1096
        _state[i] = STATE_LOWER;
1095 1097
      }
1096 1098

	
1097 1099
      // Set data for the artificial root node
1098 1100
      _root = _node_num;
1099 1101
      _parent[_root] = -1;
1100 1102
      _pred[_root] = -1;
1101 1103
      _thread[_root] = 0;
1102 1104
      _rev_thread[0] = _root;
1103 1105
      _succ_num[_root] = _node_num + 1;
1104 1106
      _last_succ[_root] = _root - 1;
1105 1107
      _supply[_root] = -_sum_supply;
1106 1108
      _pi[_root] = 0;
1107 1109

	
1108 1110
      // Add artificial arcs and initialize the spanning tree data structure
1109 1111
      if (_sum_supply == 0) {
1110 1112
        // EQ supply constraints
1111 1113
        _search_arc_num = _arc_num;
1112 1114
        _all_arc_num = _arc_num + _node_num;
1113 1115
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1114 1116
          _parent[u] = _root;
1115 1117
          _pred[u] = e;
1116 1118
          _thread[u] = u + 1;
1117 1119
          _rev_thread[u + 1] = u;
1118 1120
          _succ_num[u] = 1;
1119 1121
          _last_succ[u] = u;
1120 1122
          _cap[e] = INF;
1121 1123
          _state[e] = STATE_TREE;
1122 1124
          if (_supply[u] >= 0) {
1123 1125
            _pred_dir[u] = DIR_UP;
1124 1126
            _pi[u] = 0;
1125 1127
            _source[e] = u;
1126 1128
            _target[e] = _root;
1127 1129
            _flow[e] = _supply[u];
1128 1130
            _cost[e] = 0;
1129 1131
          } else {
1130 1132
            _pred_dir[u] = DIR_DOWN;
1131 1133
            _pi[u] = ART_COST;
1132 1134
            _source[e] = _root;
1133 1135
            _target[e] = u;
1134 1136
            _flow[e] = -_supply[u];
1135 1137
            _cost[e] = ART_COST;
1136 1138
          }
1137 1139
        }
0 comments (0 inline)