gravatar
alpar (Alpar Juttner)
alpar@cs.elte.hu
Merge
0 16 1
merge default
2 files changed:
↑ Collapse diff ↑
Ignore white space 6 line context
1
%!PS-Adobe-2.0 EPSF-2.0
2
%%Creator: LEMON, graphToEps()
3
%%CreationDate: Fri Oct 19 18:32:32 2007
4
%%BoundingBox: 0 0 596 842
5
%%DocumentPaperSizes: a4
6
%%EndComments
7
/lb { setlinewidth setrgbcolor newpath moveto
8
      4 2 roll 1 index 1 index curveto stroke } bind def
9
/l { setlinewidth setrgbcolor newpath moveto lineto stroke } bind def
10
/c { newpath dup 3 index add 2 index moveto 0 360 arc closepath } bind def
11
/sq { newpath 2 index 1 index add 2 index 2 index add moveto
12
      2 index 1 index sub 2 index 2 index add lineto
13
      2 index 1 index sub 2 index 2 index sub lineto
14
      2 index 1 index add 2 index 2 index sub lineto
15
      closepath pop pop pop} bind def
16
/di { newpath 2 index 1 index add 2 index moveto
17
      2 index             2 index 2 index add lineto
18
      2 index 1 index sub 2 index             lineto
19
      2 index             2 index 2 index sub lineto
20
      closepath pop pop pop} bind def
21
/nc { 0 0 0 setrgbcolor 5 index 5 index 5 index c fill
22
     setrgbcolor 1.1 div c fill
23
   } bind def
24
/nsq { 0 0 0 setrgbcolor 5 index 5 index 5 index sq fill
25
     setrgbcolor 1.1 div sq fill
26
   } bind def
27
/ndi { 0 0 0 setrgbcolor 5 index 5 index 5 index di fill
28
     setrgbcolor 1.1 div di fill
29
   } bind def
30
/nfemale { 0 0 0 setrgbcolor 3 index 0.0909091 1.5 mul mul setlinewidth
31
  newpath 5 index 5 index moveto 5 index 5 index 5 index 3.01 mul sub
32
  lineto 5 index 4 index .7 mul sub 5 index 5 index 2.2 mul sub moveto
33
  5 index 4 index .7 mul add 5 index 5 index 2.2 mul sub lineto stroke
34
  5 index 5 index 5 index c fill
35
  setrgbcolor 1.1 div c fill
36
  } bind def
37
/nmale {
38
  0 0 0 setrgbcolor 3 index 0.0909091 1.5 mul mul setlinewidth
39
  newpath 5 index 5 index moveto
40
  5 index 4 index 1 mul 1.5 mul add
41
  5 index 5 index 3 sqrt 1.5 mul mul add
42
  1 index 1 index lineto
43
  1 index 1 index 7 index sub moveto
44
  1 index 1 index lineto
45
  exch 5 index 3 sqrt .5 mul mul sub exch 5 index .5 mul sub lineto
46
  stroke
47
  5 index 5 index 5 index c fill
48
  setrgbcolor 1.1 div c fill
49
  } bind def
50
/arrl 1 def
51
/arrw 0.3 def
52
/lrl { 2 index mul exch 2 index mul exch rlineto pop} bind def
53
/arr { setrgbcolor /y1 exch def /x1 exch def /dy exch def /dx exch def
54
       /w exch def /len exch def
55
       newpath x1 dy w 2 div mul add y1 dx w 2 div mul sub moveto
56
       len w sub arrl sub dx dy lrl
57
       arrw dy dx neg lrl
58
       dx arrl w add mul dy w 2 div arrw add mul sub
59
       dy arrl w add mul dx w 2 div arrw add mul add rlineto
60
       dx arrl w add mul neg dy w 2 div arrw add mul sub
61
       dy arrl w add mul neg dx w 2 div arrw add mul add rlineto
62
       arrw dy dx neg lrl
63
       len w sub arrl sub neg dx dy lrl
64
       closepath fill } bind def
65
/cshow { 2 index 2 index moveto dup stringwidth pop
66
         neg 2 div fosi .35 mul neg rmoveto show pop pop} def
67

	
68
gsave
69
15 138.307 translate
70
12.7843 dup scale
71
90 rotate
72
0.608112 -43.6081 translate
73
%Edges:
74
gsave
75
9 31 9.5 30.5 10 30 0 0 0 0.091217 lb
76
9 31 5.5 34.5 2 38 0 0 0 0.091217 lb
77
9 31 25.5 16 42 1 0 0 0 0.091217 lb
78
3 40 23 20.5 43 1 0 0 0 0.091217 lb
79
3 40 22.5 20.5 42 1 0 0 0 0.091217 lb
80
3 40 2.5 40.5 2 41 0 0 0 0.091217 lb
81
13 25 10.5 24.5 8 24 0 0 0 0.091217 lb
82
13 25 12 27 11 29 0 0 0 0.091217 lb
83
3 4 2.5 3 2 2 0 0 0 0.091217 lb
84
3 4 4.5 3 6 2 0 0 0 0.091217 lb
85
6 25 7 24.5 8 24 0 0 0 0.091217 lb
86
6 25 6 24.5 6 24 0 0 0 0.091217 lb
87
34 2 33.5 2 33 2 0 0 0 0.091217 lb
88
34 2 35 2 36 2 0 0 0 0.091217 lb
89
6 8 16 9 26 10 0 0 0 0.091217 lb
90
6 8 6 10.5 6 13 0 0 0 0.091217 lb
91
6 8 6 7.5 6 7 0 0 0 0.091217 lb
92
26 10 27.5 8.5 29 7 0 0 0 0.091217 lb
93
26 10 27.5 9 29 8 0 0 0 0.091217 lb
94
10 30 10.5 29.5 11 29 0 0 0 0.091217 lb
95
8 24 7 23.5 6 23 0 0 0 0.091217 lb
96
8 24 8 24.5 8 25 0 0 0 0.091217 lb
97
33 2 32.5 2 32 2 0 0 0 0.091217 lb
98
29 7 17.5 7 6 7 0 0 0 0.091217 lb
99
2 2 1.5 22 1 42 0 0 0 0.091217 lb
100
2 2 3.5 2 5 2 0 0 0 0.091217 lb
101
21 15 13.5 14.5 6 14 0 0 0 0.091217 lb
102
21 15 21 15.5 21 16 0 0 0 0.091217 lb
103
1 42 0.5 42.5 0 43 0 0 0 0.091217 lb
104
1 42 1.5 41.5 2 41 0 0 0 0.091217 lb
105
6 15 6 15.5 6 16 0 0 0 0.091217 lb
106
6 15 6 14.5 6 14 0 0 0 0.091217 lb
107
43 1 22 0.5 1 0 0 0 0 0.091217 lb
108
31 2 18.5 2 6 2 0 0 0 0.091217 lb
109
31 2 31.5 2 32 2 0 0 0 0.091217 lb
110
6 24 6 23.5 6 23 0 0 0 0.091217 lb
111
6 16 6 16.5 6 17 0 0 0 0.091217 lb
112
6 23 6 20 6 17 0 0 0 0.091217 lb
113
6 2 5.5 2 5 2 0 0 0 0.091217 lb
114
6 2 6 4.5 6 7 0 0 0 0.091217 lb
115
0 43 0.5 21.5 1 0 0 0 0 0.091217 lb
116
1 1 19.5 1.5 38 2 0 0 0 0.091217 lb
117
1 1 1 0.5 1 0 0 0 0 0.091217 lb
118
2 38 5.5 31.5 9 25 0 0 0 0.091217 lb
119
25 13 15.5 13 6 13 0 0 0 0.091217 lb
120
25 13 15.5 13.5 6 14 0 0 0 0.091217 lb
121
8 25 8.5 25 9 25 0 0 0 0.091217 lb
122
11 29 24.5 15.5 38 2 0 0 0 0.091217 lb
123
6 17 11.5 18 17 19 0 0 0 0.091217 lb
124
16 23 26.5 12.5 37 2 0 0 0 0.091217 lb
125
16 23 18.5 19.5 21 16 0 0 0 0.091217 lb
126
36 2 36.5 2 37 2 0 0 0 0.091217 lb
127
36 2 32.5 5 29 8 0 0 0 0.091217 lb
128
6 13 6 13.5 6 14 0 0 0 0.091217 lb
129
37 2 37.5 2 38 2 0 0 0 0.091217 lb
130
21 16 19 17.5 17 19 0 0 0 0.091217 lb
131
grestore
132
%Nodes:
133
gsave
134
29 8 0.304556 1 1 1 nc
135
2 41 0.304556 1 1 1 nc
136
6 7 0.304556 1 1 1 nc
137
5 2 0.304556 1 1 1 nc
138
17 19 0.304556 1 1 1 nc
139
21 16 0.304556 1 1 1 nc
140
1 0 0.304556 1 1 1 nc
141
9 25 0.304556 1 1 1 nc
142
6 14 0.304556 1 1 1 nc
143
42 1 0.304556 1 1 1 nc
144
38 2 0.304556 1 1 1 nc
145
37 2 0.304556 1 1 1 nc
146
6 13 0.304556 1 1 1 nc
147
36 2 0.304556 1 1 1 nc
148
16 23 0.304556 1 1 1 nc
149
6 17 0.304556 1 1 1 nc
150
11 29 0.304556 1 1 1 nc
151
8 25 0.304556 1 1 1 nc
152
32 2 0.304556 1 1 1 nc
153
25 13 0.304556 1 1 1 nc
154
2 38 0.304556 1 1 1 nc
155
1 1 0.304556 1 1 1 nc
156
0 43 0.304556 1 1 1 nc
157
6 2 0.304556 1 1 1 nc
158
6 23 0.304556 1 1 1 nc
159
6 16 0.304556 1 1 1 nc
160
6 24 0.304556 1 1 1 nc
161
31 2 0.304556 1 1 1 nc
162
43 1 0.304556 1 1 1 nc
163
6 15 0.304556 1 1 1 nc
164
1 42 0.304556 1 1 1 nc
165
21 15 0.304556 1 1 1 nc
166
2 2 0.304556 1 1 1 nc
167
29 7 0.304556 1 1 1 nc
168
33 2 0.304556 1 1 1 nc
169
8 24 0.304556 1 1 1 nc
170
10 30 0.304556 1 1 1 nc
171
26 10 0.304556 1 1 1 nc
172
6 8 0.304556 1 1 1 nc
173
34 2 0.304556 1 1 1 nc
174
6 25 0.304556 1 1 1 nc
175
3 4 0.304556 1 1 1 nc
176
13 25 0.304556 1 1 1 nc
177
3 40 0.304556 1 1 1 nc
178
9 31 0.304556 1 1 1 nc
179
grestore
180
grestore
181
showpage
Ignore white space 6 line context
1 1
Installation Instructions
2 2
=========================
3 3

	
4 4
Since you are reading this I assume you already obtained one of the release
5 5
tarballs and successfully extracted it. The latest version of LEMON is
6 6
available at our web page (http://lemon.cs.elte.hu/).
7 7

	
8 8
LEMON provides two different build environments, one is based on "autotool",
9 9
while the other is based on "cmake". This file contains instructions only for
10 10
the former one, which is the recommended build environment on Linux, Mac OSX
11 11
and other unices or if you use Cygwin on Windows. For cmake installation
12 12
instructions visit http://lemon.cs.elte.hu.
13 13

	
14 14
In order to install LEMON from the extracted source tarball you have to
15 15
issue the following commands:
16 16

	
17 17
   1. `cd lemon-x.y.z'
18 18

	
19 19
      This command changes to the directory which was created when you
20 20
      extracted the sources. The x.y.z part is a version number.
21 21

	
22 22
   2. `./configure'
23 23

	
24 24
      This command runs the configure shell script, which does some checks and
25 25
      creates the makefiles.
26 26

	
27 27
   3. `make'
28 28

	
29 29
      This command compiles the non-template part of LEMON into libemon.a
30 30
      file. It also compiles the programs in the tools subdirectory by
31 31
      default.
32 32

	
33 33
   4. `make check'
34 34

	
35 35
      This step is optional, but recommended. It runs the test programs that
36 36
      we developed for LEMON to check whether the library works properly on
37 37
      your platform.
38 38

	
39 39
   5. `make install'
40 40

	
41 41
      This command installs LEMON under /usr/local (you will need root
42 42
      privileges to be able to do that). If you want to install it to some
43 43
      other location, then pass the --prefix=DIRECTORY flag to configure in
44 44
      step 2. For example: `./configure --prefix=/home/username/lemon'.
45 45

	
46 46
   6. `make install-html'
47 47

	
48 48
      This command installs the documentation under share/doc/lemon/docs. The
49 49
      generated documentation is included in the tarball. If you want to
50 50
      generate it yourself, then run `make html'. Note that for this you need
51 51
      to have the following programs installed: Doxygen, Graphviz, Ghostscript,
52 52
      Latex.
53 53

	
54 54

	
55 55
Configure Options and Variables
56 56
===============================
57 57

	
58 58
In step 2 you can customize the actions of configure by setting variables
59 59
and passing options to it. This can be done like this:
60 60
`./configure [OPTION]... [VARIABLE=VALUE]...'
61 61

	
62 62
Below you will find some useful variables and options (see `./configure --help'
63 63
for more):
64 64

	
65 65
CXX='comp'
66 66

	
67 67
  Change the C++ compiler to 'comp'.
68 68

	
69 69
CXXFLAGS='flags'
70 70

	
71 71
  Pass the 'flags' to the compiler. For example CXXFLAGS='-O3 -march=pentium-m'
72 72
  turns on generation of aggressively optimized Pentium-M specific code.
73 73

	
74 74
--prefix=PREFIX
75 75

	
76 76
  Set the installation prefix to PREFIX. By default it is /usr/local.
77 77

	
78 78
--enable-tools
79 79

	
80 80
   Build the programs in the tools subdirectory (default).
81 81

	
82 82
--disable-tools
83 83

	
84 84
   Do not build the programs in the tools subdirectory.
85 85

	
86 86
--with-glpk[=PREFIX]
87 87

	
88 88
   Enable GLPK support (default). You should specify the prefix too if
89 89
   you installed GLPK to some non-standard location (e.g. your home
90 90
   directory). If it is not found, GLPK support will be disabled.
91 91

	
92 92
--with-glpk-includedir=DIR
93 93

	
94 94
   The directory where the GLPK header files are located. This is only
95 95
   useful when the GLPK headers and libraries are not under the same
96 96
   prefix (which is unlikely).
97 97

	
98 98
--with-glpk-libdir=DIR
99 99

	
100 100
   The directory where the GLPK libraries are located. This is only
101 101
   useful when the GLPK headers and libraries are not under the same
102 102
   prefix (which is unlikely).
103 103

	
104 104
--without-glpk
105 105

	
106 106
   Disable GLPK support.
107 107

	
108 108
--with-cplex[=PREFIX]
109 109

	
110 110
   Enable CPLEX support (default). You should specify the prefix too
111 111
   if you installed CPLEX to some non-standard location
112 112
   (e.g. /opt/ilog/cplex75). If it is not found, CPLEX support will be
113 113
   disabled.
114 114

	
115 115
--with-cplex-includedir=DIR
116 116

	
117 117
   The directory where the CPLEX header files are located. This is
118 118
   only useful when the CPLEX headers and libraries are not under the
119 119
   same prefix (e.g.  /usr/local/cplex/cplex75/include).
120 120

	
121 121
--with-cplex-libdir=DIR
122 122

	
123 123
   The directory where the CPLEX libraries are located. This is only
124 124
   useful when the CPLEX headers and libraries are not under the same
125 125
   prefix (e.g.
126 126
   /usr/local/cplex/cplex75/lib/i86_linux2_glibc2.2_gcc3.0/static_pic_mt).
127 127

	
128 128
--without-cplex
129 129

	
130 130
   Disable CPLEX support.
131 131

	
132 132
--with-soplex[=PREFIX]
133 133

	
134 134
   Enable SoPlex support (default). You should specify the prefix too if
135 135
   you installed SoPlex to some non-standard location (e.g. your home
136 136
   directory). If it is not found, SoPlex support will be disabled.
137 137

	
138 138
--with-soplex-includedir=DIR
139 139

	
140 140
   The directory where the SoPlex header files are located. This is only
141 141
   useful when the SoPlex headers and libraries are not under the same
142 142
   prefix (which is unlikely).
143 143

	
144 144
--with-soplex-libdir=DIR
145 145

	
146 146
   The directory where the SoPlex libraries are located. This is only
147 147
   useful when the SoPlex headers and libraries are not under the same
148 148
   prefix (which is unlikely).
149 149

	
150 150
--without-soplex
151 151

	
152 152
   Disable SoPlex support.
153 153

	
154 154
--with-coin[=PREFIX]
155 155

	
156 156
   Enable support for COIN-OR solvers (CLP and CBC). You should
157 157
   specify the prefix too. (by default, COIN-OR tools install
158 158
   themselves to the source code directory). This command enables the
159 159
   solvers that are actually found.
160 160

	
161 161
--with-coin-includedir=DIR
162 162

	
163 163
   The directory where the COIN-OR header files are located. This is
164 164
   only useful when the COIN-OR headers and libraries are not under
165 165
   the same prefix (which is unlikely).
166 166

	
167 167
--with-coin-libdir=DIR
168 168

	
169 169
   The directory where the COIN-OR libraries are located. This is only
170 170
   useful when the COIN-OR headers and libraries are not under the
171 171
   same prefix (which is unlikely).
172 172

	
173 173
--without-coin
174 174

	
175 175
   Disable COIN-OR support.
176

	
177

	
178
Makefile Variables
179
==================
180

	
181
Some Makefile variables are reserved by the GNU Coding Standards for
182
the use of the "user" - the person building the package. For instance,
183
CXX and CXXFLAGS are such variables, and have the same meaning as
184
explained in the previous section. These variables can be set on the
185
command line when invoking `make' like this:
186
`make [VARIABLE=VALUE]...'
187

	
188
WARNINGCXXFLAGS is a non-standard Makefile variable introduced by us
189
to hold several compiler flags related to warnings. Its default value
190
can be overridden when invoking `make'. For example to disable all
191
warning flags use `make WARNINGCXXFLAGS='.
192

	
193
In order to turn off a single flag from the default set of warning
194
flags, you can use the CXXFLAGS variable, since this is passed after
195
WARNINGCXXFLAGS. For example to turn off `-Wold-style-cast' (which is
196
used by default when g++ is detected) you can use
197
`make CXXFLAGS="-g -O2 -Wno-old-style-cast"'.
Ignore white space 6 line context
1 1
SET(PACKAGE_NAME ${PROJECT_NAME})
2 2
SET(PACKAGE_VERSION ${PROJECT_VERSION})
3 3
SET(abs_top_srcdir ${PROJECT_SOURCE_DIR})
4 4
SET(abs_top_builddir ${PROJECT_BINARY_DIR})
5 5

	
6 6
CONFIGURE_FILE(
7 7
  ${PROJECT_SOURCE_DIR}/doc/Doxyfile.in
8 8
  ${PROJECT_BINARY_DIR}/doc/Doxyfile
9 9
  @ONLY
10 10
)
11 11

	
12 12
IF(DOXYGEN_EXECUTABLE AND PYTHONINTERP_FOUND AND GHOSTSCRIPT_EXECUTABLE)
13 13
  FILE(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/html/)
14 14
  SET(GHOSTSCRIPT_OPTIONS -dNOPAUSE -dBATCH -q -dEPSCrop -dTextAlphaBits=4 -dGraphicsAlphaBits=4 -sDEVICE=pngalpha)
15 15
  ADD_CUSTOM_TARGET(html
16 16
    COMMAND ${CMAKE_COMMAND} -E remove_directory gen-images
17 17
    COMMAND ${CMAKE_COMMAND} -E make_directory gen-images
18 18
    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/bipartite_matching.png ${CMAKE_CURRENT_SOURCE_DIR}/images/bipartite_matching.eps
19 19
    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/bipartite_partitions.png ${CMAKE_CURRENT_SOURCE_DIR}/images/bipartite_partitions.eps
20 20
    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/connected_components.png ${CMAKE_CURRENT_SOURCE_DIR}/images/connected_components.eps
21 21
    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/edge_biconnected_components.png ${CMAKE_CURRENT_SOURCE_DIR}/images/edge_biconnected_components.eps
22 22
    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/grid_graph.png ${CMAKE_CURRENT_SOURCE_DIR}/images/grid_graph.eps
23 23
    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/node_biconnected_components.png ${CMAKE_CURRENT_SOURCE_DIR}/images/node_biconnected_components.eps
24 24
    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/nodeshape_0.png ${CMAKE_CURRENT_SOURCE_DIR}/images/nodeshape_0.eps
25 25
    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/nodeshape_1.png ${CMAKE_CURRENT_SOURCE_DIR}/images/nodeshape_1.eps
26 26
    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/nodeshape_2.png ${CMAKE_CURRENT_SOURCE_DIR}/images/nodeshape_2.eps
27 27
    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/nodeshape_3.png ${CMAKE_CURRENT_SOURCE_DIR}/images/nodeshape_3.eps
28 28
    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/nodeshape_4.png ${CMAKE_CURRENT_SOURCE_DIR}/images/nodeshape_4.eps
29
    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/planar.png ${CMAKE_CURRENT_SOURCE_DIR}/images/planar.eps
29 30
    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/strongly_connected_components.png ${CMAKE_CURRENT_SOURCE_DIR}/images/strongly_connected_components.eps
30 31
    COMMAND ${CMAKE_COMMAND} -E remove_directory html
31 32
    COMMAND ${PYTHON_EXECUTABLE} ${PROJECT_SOURCE_DIR}/scripts/bib2dox.py ${CMAKE_CURRENT_SOURCE_DIR}/references.bib >references.dox
32 33
    COMMAND ${DOXYGEN_EXECUTABLE} Doxyfile
33 34
    WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
34 35
  )
35 36

	
36 37
  SET_TARGET_PROPERTIES(html PROPERTIES PROJECT_LABEL BUILD_DOC)
37 38

	
38 39
  IF(UNIX)
39 40
    INSTALL(
40 41
      DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/html/
41 42
      DESTINATION share/doc/lemon/html
42 43
      COMPONENT html_documentation
43 44
    )
44 45
  ELSEIF(WIN32)
45 46
    INSTALL(
46 47
      DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/html/
47 48
      DESTINATION doc
48 49
      COMPONENT html_documentation
49 50
    )
50 51
  ENDIF()
51 52

	
52 53
ENDIF()
Ignore white space 6 line context
1 1
EXTRA_DIST += \
2 2
	doc/Doxyfile.in \
3 3
	doc/DoxygenLayout.xml \
4 4
	doc/coding_style.dox \
5 5
	doc/dirs.dox \
6 6
	doc/groups.dox \
7 7
	doc/lgf.dox \
8 8
	doc/license.dox \
9 9
	doc/mainpage.dox \
10 10
	doc/migration.dox \
11 11
	doc/min_cost_flow.dox \
12 12
	doc/named-param.dox \
13 13
	doc/namespaces.dox \
14 14
	doc/html \
15 15
	doc/CMakeLists.txt
16 16

	
17 17
DOC_EPS_IMAGES18 = \
18 18
	grid_graph.eps \
19 19
	nodeshape_0.eps \
20 20
	nodeshape_1.eps \
21 21
	nodeshape_2.eps \
22 22
	nodeshape_3.eps \
23 23
	nodeshape_4.eps
24 24

	
25 25
DOC_EPS_IMAGES27 = \
26 26
	bipartite_matching.eps \
27 27
	bipartite_partitions.eps \
28 28
	connected_components.eps \
29 29
	edge_biconnected_components.eps \
30 30
	node_biconnected_components.eps \
31
	planar.eps \
31 32
	strongly_connected_components.eps
32 33

	
33 34
DOC_EPS_IMAGES = \
34 35
	$(DOC_EPS_IMAGES18) \
35 36
	$(DOC_EPS_IMAGES27)
36 37

	
37 38
DOC_PNG_IMAGES = \
38 39
	$(DOC_EPS_IMAGES:%.eps=doc/gen-images/%.png)
39 40

	
40 41
EXTRA_DIST += $(DOC_EPS_IMAGES:%=doc/images/%)
41 42

	
42 43
doc/html:
43 44
	$(MAKE) $(AM_MAKEFLAGS) html
44 45

	
45 46
GS_COMMAND=gs -dNOPAUSE -dBATCH -q -dEPSCrop -dTextAlphaBits=4 -dGraphicsAlphaBits=4
46 47

	
47 48
$(DOC_EPS_IMAGES18:%.eps=doc/gen-images/%.png): doc/gen-images/%.png: doc/images/%.eps
48 49
	-mkdir doc/gen-images
49 50
	if test ${gs_found} = yes; then \
50 51
	  $(GS_COMMAND) -sDEVICE=pngalpha -r18 -sOutputFile=$@ $<; \
51 52
	else \
52 53
	  echo; \
53 54
	  echo "Ghostscript not found."; \
54 55
	  echo; \
55 56
	  exit 1; \
56 57
	fi
57 58

	
58 59
$(DOC_EPS_IMAGES27:%.eps=doc/gen-images/%.png): doc/gen-images/%.png: doc/images/%.eps
59 60
	-mkdir doc/gen-images
60 61
	if test ${gs_found} = yes; then \
61 62
	  $(GS_COMMAND) -sDEVICE=pngalpha -r27 -sOutputFile=$@ $<; \
62 63
	else \
63 64
	  echo; \
64 65
	  echo "Ghostscript not found."; \
65 66
	  echo; \
66 67
	  exit 1; \
67 68
	fi
68 69

	
69 70
references.dox: doc/references.bib
70 71
	if test ${python_found} = yes; then \
71 72
	  cd doc; \
72 73
	  python @abs_top_srcdir@/scripts/bib2dox.py @abs_top_builddir@/$< >$@; \
73 74
	  cd ..; \
74 75
	else \
75 76
	  echo; \
76 77
	  echo "Python not found."; \
77 78
	  echo; \
78 79
	  exit 1; \
79 80
	fi
80 81

	
81 82
html-local: $(DOC_PNG_IMAGES) references.dox
82 83
	if test ${doxygen_found} = yes; then \
83 84
	  cd doc; \
84 85
	  doxygen Doxyfile; \
85 86
	  cd ..; \
86 87
	else \
87 88
	  echo; \
88 89
	  echo "Doxygen not found."; \
89 90
	  echo; \
90 91
	  exit 1; \
91 92
	fi
92 93

	
93 94
clean-local:
94 95
	-rm -rf doc/html
95 96
	-rm -f doc/doxygen.log
96 97
	-rm -f $(DOC_PNG_IMAGES)
97 98
	-rm -rf doc/gen-images
98 99

	
99 100
update-external-tags:
100 101
	wget -O doc/libstdc++.tag.tmp http://gcc.gnu.org/onlinedocs/libstdc++/latest-doxygen/libstdc++.tag && \
101 102
	mv doc/libstdc++.tag.tmp doc/libstdc++.tag || \
102 103
	rm doc/libstdc++.tag.tmp
103 104

	
104 105
install-html-local: doc/html
105 106
	@$(NORMAL_INSTALL)
106 107
	$(mkinstalldirs) $(DESTDIR)$(htmldir)/html
107 108
	for p in doc/html/*.{html,css,png,map,gif,tag} ; do \
108 109
	  f="`echo $$p | sed -e 's|^.*/||'`"; \
109 110
	  echo " $(INSTALL_DATA) $$p $(DESTDIR)$(htmldir)/html/$$f"; \
110 111
	  $(INSTALL_DATA) $$p $(DESTDIR)$(htmldir)/html/$$f; \
111 112
	done
112 113

	
113 114
uninstall-local:
114 115
	@$(NORMAL_UNINSTALL)
115 116
	for p in doc/html/*.{html,css,png,map,gif,tag} ; do \
116 117
	  f="`echo $$p | sed -e 's|^.*/||'`"; \
117 118
	  echo " rm -f $(DESTDIR)$(htmldir)/html/$$f"; \
118 119
	  rm -f $(DESTDIR)$(htmldir)/html/$$f; \
119 120
	done
120 121

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

	
19 19
#ifndef LEMON_BELLMAN_FORD_H
20 20
#define LEMON_BELLMAN_FORD_H
21 21

	
22 22
/// \ingroup shortest_path
23 23
/// \file
24 24
/// \brief Bellman-Ford algorithm.
25 25

	
26 26
#include <lemon/list_graph.h>
27 27
#include <lemon/bits/path_dump.h>
28 28
#include <lemon/core.h>
29 29
#include <lemon/error.h>
30 30
#include <lemon/maps.h>
31 31
#include <lemon/path.h>
32 32

	
33 33
#include <limits>
34 34

	
35 35
namespace lemon {
36 36

	
37 37
  /// \brief Default OperationTraits for the BellmanFord algorithm class.
38 38
  ///  
39 39
  /// This operation traits class defines all computational operations
40 40
  /// and constants that are used in the Bellman-Ford algorithm.
41 41
  /// The default implementation is based on the \c numeric_limits class.
42 42
  /// If the numeric type does not have infinity value, then the maximum
43 43
  /// value is used as extremal infinity value.
44 44
  template <
45 45
    typename V, 
46 46
    bool has_inf = std::numeric_limits<V>::has_infinity>
47 47
  struct BellmanFordDefaultOperationTraits {
48 48
    /// \e
49 49
    typedef V Value;
50 50
    /// \brief Gives back the zero value of the type.
51 51
    static Value zero() {
52 52
      return static_cast<Value>(0);
53 53
    }
54 54
    /// \brief Gives back the positive infinity value of the type.
55 55
    static Value infinity() {
56 56
      return std::numeric_limits<Value>::infinity();
57 57
    }
58 58
    /// \brief Gives back the sum of the given two elements.
59 59
    static Value plus(const Value& left, const Value& right) {
60 60
      return left + right;
61 61
    }
62 62
    /// \brief Gives back \c true only if the first value is less than
63 63
    /// the second.
64 64
    static bool less(const Value& left, const Value& right) {
65 65
      return left < right;
66 66
    }
67 67
  };
68 68

	
69 69
  template <typename V>
70 70
  struct BellmanFordDefaultOperationTraits<V, false> {
71 71
    typedef V Value;
72 72
    static Value zero() {
73 73
      return static_cast<Value>(0);
74 74
    }
75 75
    static Value infinity() {
76 76
      return std::numeric_limits<Value>::max();
77 77
    }
78 78
    static Value plus(const Value& left, const Value& right) {
79 79
      if (left == infinity() || right == infinity()) return infinity();
80 80
      return left + right;
81 81
    }
82 82
    static bool less(const Value& left, const Value& right) {
83 83
      return left < right;
84 84
    }
85 85
  };
86 86
  
87 87
  /// \brief Default traits class of BellmanFord class.
88 88
  ///
89 89
  /// Default traits class of BellmanFord class.
90 90
  /// \param GR The type of the digraph.
91 91
  /// \param LEN The type of the length map.
92 92
  template<typename GR, typename LEN>
93 93
  struct BellmanFordDefaultTraits {
94 94
    /// The type of the digraph the algorithm runs on. 
95 95
    typedef GR Digraph;
96 96

	
97 97
    /// \brief The type of the map that stores the arc lengths.
98 98
    ///
99 99
    /// The type of the map that stores the arc lengths.
100 100
    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
101 101
    typedef LEN LengthMap;
102 102

	
103 103
    /// The type of the arc lengths.
104 104
    typedef typename LEN::Value Value;
105 105

	
106 106
    /// \brief Operation traits for Bellman-Ford algorithm.
107 107
    ///
108 108
    /// It defines the used operations and the infinity value for the
109 109
    /// given \c Value type.
110 110
    /// \see BellmanFordDefaultOperationTraits
111 111
    typedef BellmanFordDefaultOperationTraits<Value> OperationTraits;
112 112
 
113 113
    /// \brief The type of the map that stores the last arcs of the 
114 114
    /// shortest paths.
115 115
    /// 
116 116
    /// The type of the map that stores the last
117 117
    /// arcs of the shortest paths.
118 118
    /// It must conform to the \ref concepts::WriteMap "WriteMap" concept.
119 119
    typedef typename GR::template NodeMap<typename GR::Arc> PredMap;
120 120

	
121 121
    /// \brief Instantiates a \c PredMap.
122 122
    /// 
123 123
    /// This function instantiates a \ref PredMap. 
124 124
    /// \param g is the digraph to which we would like to define the
125 125
    /// \ref PredMap.
126 126
    static PredMap *createPredMap(const GR& g) {
127 127
      return new PredMap(g);
128 128
    }
129 129

	
130 130
    /// \brief The type of the map that stores the distances of the nodes.
131 131
    ///
132 132
    /// The type of the map that stores the distances of the nodes.
133 133
    /// It must conform to the \ref concepts::WriteMap "WriteMap" concept.
134 134
    typedef typename GR::template NodeMap<typename LEN::Value> DistMap;
135 135

	
136 136
    /// \brief Instantiates a \c DistMap.
137 137
    ///
138 138
    /// This function instantiates a \ref DistMap. 
139 139
    /// \param g is the digraph to which we would like to define the 
140 140
    /// \ref DistMap.
141 141
    static DistMap *createDistMap(const GR& g) {
142 142
      return new DistMap(g);
143 143
    }
144 144

	
145 145
  };
146 146
  
147 147
  /// \brief %BellmanFord algorithm class.
148 148
  ///
149 149
  /// \ingroup shortest_path
150 150
  /// This class provides an efficient implementation of the Bellman-Ford 
151 151
  /// algorithm. The maximum time complexity of the algorithm is
152 152
  /// <tt>O(ne)</tt>.
153 153
  ///
154 154
  /// The Bellman-Ford algorithm solves the single-source shortest path
155 155
  /// problem when the arcs can have negative lengths, but the digraph
156 156
  /// should not contain directed cycles with negative total length.
157 157
  /// If all arc costs are non-negative, consider to use the Dijkstra
158 158
  /// algorithm instead, since it is more efficient.
159 159
  ///
160 160
  /// The arc lengths are passed to the algorithm using a
161 161
  /// \ref concepts::ReadMap "ReadMap", so it is easy to change it to any 
162 162
  /// kind of length. The type of the length values is determined by the
163 163
  /// \ref concepts::ReadMap::Value "Value" type of the length map.
164 164
  ///
165 165
  /// There is also a \ref bellmanFord() "function-type interface" for the
166 166
  /// Bellman-Ford algorithm, which is convenient in the simplier cases and
167 167
  /// it can be used easier.
168 168
  ///
169 169
  /// \tparam GR The type of the digraph the algorithm runs on.
170 170
  /// The default type is \ref ListDigraph.
171 171
  /// \tparam LEN A \ref concepts::ReadMap "readable" arc map that specifies
172 172
  /// the lengths of the arcs. The default map type is
173 173
  /// \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
174
  /// \tparam TR The traits class that defines various types used by the
175
  /// algorithm. By default, it is \ref BellmanFordDefaultTraits
176
  /// "BellmanFordDefaultTraits<GR, LEN>".
177
  /// In most cases, this parameter should not be set directly,
178
  /// consider to use the named template parameters instead.
174 179
#ifdef DOXYGEN
175 180
  template <typename GR, typename LEN, typename TR>
176 181
#else
177 182
  template <typename GR=ListDigraph,
178 183
            typename LEN=typename GR::template ArcMap<int>,
179 184
            typename TR=BellmanFordDefaultTraits<GR,LEN> >
180 185
#endif
181 186
  class BellmanFord {
182 187
  public:
183 188

	
184 189
    ///The type of the underlying digraph.
185 190
    typedef typename TR::Digraph Digraph;
186 191
    
187 192
    /// \brief The type of the arc lengths.
188 193
    typedef typename TR::LengthMap::Value Value;
189 194
    /// \brief The type of the map that stores the arc lengths.
190 195
    typedef typename TR::LengthMap LengthMap;
191 196
    /// \brief The type of the map that stores the last
192 197
    /// arcs of the shortest paths.
193 198
    typedef typename TR::PredMap PredMap;
194 199
    /// \brief The type of the map that stores the distances of the nodes.
195 200
    typedef typename TR::DistMap DistMap;
196 201
    /// The type of the paths.
197 202
    typedef PredMapPath<Digraph, PredMap> Path;
198 203
    ///\brief The \ref BellmanFordDefaultOperationTraits
199 204
    /// "operation traits class" of the algorithm.
200 205
    typedef typename TR::OperationTraits OperationTraits;
201 206

	
202 207
    ///The \ref BellmanFordDefaultTraits "traits class" of the algorithm.
203 208
    typedef TR Traits;
204 209

	
205 210
  private:
206 211

	
207 212
    typedef typename Digraph::Node Node;
208 213
    typedef typename Digraph::NodeIt NodeIt;
209 214
    typedef typename Digraph::Arc Arc;
210 215
    typedef typename Digraph::OutArcIt OutArcIt;
211 216

	
212 217
    // Pointer to the underlying digraph.
213 218
    const Digraph *_gr;
214 219
    // Pointer to the length map
215 220
    const LengthMap *_length;
216 221
    // Pointer to the map of predecessors arcs.
217 222
    PredMap *_pred;
218 223
    // Indicates if _pred is locally allocated (true) or not.
219 224
    bool _local_pred;
220 225
    // Pointer to the map of distances.
221 226
    DistMap *_dist;
222 227
    // Indicates if _dist is locally allocated (true) or not.
223 228
    bool _local_dist;
224 229

	
225 230
    typedef typename Digraph::template NodeMap<bool> MaskMap;
226 231
    MaskMap *_mask;
227 232

	
228 233
    std::vector<Node> _process;
229 234

	
230 235
    // Creates the maps if necessary.
231 236
    void create_maps() {
232 237
      if(!_pred) {
233 238
	_local_pred = true;
234 239
	_pred = Traits::createPredMap(*_gr);
235 240
      }
236 241
      if(!_dist) {
237 242
	_local_dist = true;
238 243
	_dist = Traits::createDistMap(*_gr);
239 244
      }
240 245
      if(!_mask) {
241 246
        _mask = new MaskMap(*_gr);
242 247
      }
243 248
    }
244 249
    
245 250
  public :
246 251
 
247 252
    typedef BellmanFord Create;
248 253

	
249 254
    /// \name Named Template Parameters
250 255

	
251 256
    ///@{
252 257

	
253 258
    template <class T>
254 259
    struct SetPredMapTraits : public Traits {
255 260
      typedef T PredMap;
256 261
      static PredMap *createPredMap(const Digraph&) {
257 262
        LEMON_ASSERT(false, "PredMap is not initialized");
258 263
        return 0; // ignore warnings
259 264
      }
260 265
    };
261 266

	
262 267
    /// \brief \ref named-templ-param "Named parameter" for setting
263 268
    /// \c PredMap type.
264 269
    ///
265 270
    /// \ref named-templ-param "Named parameter" for setting
266 271
    /// \c PredMap type.
267 272
    /// It must conform to the \ref concepts::WriteMap "WriteMap" concept.
268 273
    template <class T>
269 274
    struct SetPredMap 
270 275
      : public BellmanFord< Digraph, LengthMap, SetPredMapTraits<T> > {
271 276
      typedef BellmanFord< Digraph, LengthMap, SetPredMapTraits<T> > Create;
272 277
    };
273 278
    
274 279
    template <class T>
275 280
    struct SetDistMapTraits : public Traits {
276 281
      typedef T DistMap;
277 282
      static DistMap *createDistMap(const Digraph&) {
278 283
        LEMON_ASSERT(false, "DistMap is not initialized");
279 284
        return 0; // ignore warnings
280 285
      }
281 286
    };
282 287

	
283 288
    /// \brief \ref named-templ-param "Named parameter" for setting
284 289
    /// \c DistMap type.
285 290
    ///
286 291
    /// \ref named-templ-param "Named parameter" for setting
287 292
    /// \c DistMap type.
288 293
    /// It must conform to the \ref concepts::WriteMap "WriteMap" concept.
289 294
    template <class T>
290 295
    struct SetDistMap 
291 296
      : public BellmanFord< Digraph, LengthMap, SetDistMapTraits<T> > {
292 297
      typedef BellmanFord< Digraph, LengthMap, SetDistMapTraits<T> > Create;
293 298
    };
294 299

	
295 300
    template <class T>
296 301
    struct SetOperationTraitsTraits : public Traits {
297 302
      typedef T OperationTraits;
298 303
    };
299 304
    
300 305
    /// \brief \ref named-templ-param "Named parameter" for setting 
301 306
    /// \c OperationTraits type.
302 307
    ///
303 308
    /// \ref named-templ-param "Named parameter" for setting
304 309
    /// \c OperationTraits type.
305 310
    /// For more information, see \ref BellmanFordDefaultOperationTraits.
306 311
    template <class T>
307 312
    struct SetOperationTraits
308 313
      : public BellmanFord< Digraph, LengthMap, SetOperationTraitsTraits<T> > {
309 314
      typedef BellmanFord< Digraph, LengthMap, SetOperationTraitsTraits<T> >
310 315
      Create;
311 316
    };
312 317
    
313 318
    ///@}
314 319

	
315 320
  protected:
316 321
    
317 322
    BellmanFord() {}
318 323

	
319 324
  public:      
320 325
    
321 326
    /// \brief Constructor.
322 327
    ///
323 328
    /// Constructor.
324 329
    /// \param g The digraph the algorithm runs on.
325 330
    /// \param length The length map used by the algorithm.
326 331
    BellmanFord(const Digraph& g, const LengthMap& length) :
327 332
      _gr(&g), _length(&length),
328 333
      _pred(0), _local_pred(false),
329 334
      _dist(0), _local_dist(false), _mask(0) {}
330 335
    
331 336
    ///Destructor.
332 337
    ~BellmanFord() {
333 338
      if(_local_pred) delete _pred;
334 339
      if(_local_dist) delete _dist;
335 340
      if(_mask) delete _mask;
336 341
    }
337 342

	
338 343
    /// \brief Sets the length map.
339 344
    ///
340 345
    /// Sets the length map.
341 346
    /// \return <tt>(*this)</tt>
342 347
    BellmanFord &lengthMap(const LengthMap &map) {
343 348
      _length = &map;
344 349
      return *this;
345 350
    }
346 351

	
347 352
    /// \brief Sets the map that stores the predecessor arcs.
348 353
    ///
349 354
    /// Sets the map that stores the predecessor arcs.
350 355
    /// If you don't use this function before calling \ref run()
351 356
    /// or \ref init(), an instance will be allocated automatically.
352 357
    /// The destructor deallocates this automatically allocated map,
353 358
    /// of course.
354 359
    /// \return <tt>(*this)</tt>
355 360
    BellmanFord &predMap(PredMap &map) {
356 361
      if(_local_pred) {
357 362
	delete _pred;
358 363
	_local_pred=false;
359 364
      }
360 365
      _pred = &map;
361 366
      return *this;
362 367
    }
363 368

	
364 369
    /// \brief Sets the map that stores the distances of the nodes.
365 370
    ///
366 371
    /// Sets the map that stores the distances of the nodes calculated
367 372
    /// by the algorithm.
368 373
    /// If you don't use this function before calling \ref run()
369 374
    /// or \ref init(), an instance will be allocated automatically.
370 375
    /// The destructor deallocates this automatically allocated map,
371 376
    /// of course.
372 377
    /// \return <tt>(*this)</tt>
373 378
    BellmanFord &distMap(DistMap &map) {
374 379
      if(_local_dist) {
375 380
	delete _dist;
376 381
	_local_dist=false;
377 382
      }
378 383
      _dist = &map;
379 384
      return *this;
380 385
    }
381 386

	
382 387
    /// \name Execution Control
383 388
    /// The simplest way to execute the Bellman-Ford algorithm is to use
384 389
    /// one of the member functions called \ref run().\n
385 390
    /// If you need better control on the execution, you have to call
386 391
    /// \ref init() first, then you can add several source nodes
387 392
    /// with \ref addSource(). Finally the actual path computation can be
388 393
    /// performed with \ref start(), \ref checkedStart() or
389 394
    /// \ref limitedStart().
390 395

	
391 396
    ///@{
392 397

	
393 398
    /// \brief Initializes the internal data structures.
394 399
    /// 
395 400
    /// Initializes the internal data structures. The optional parameter
396 401
    /// is the initial distance of each node.
397 402
    void init(const Value value = OperationTraits::infinity()) {
398 403
      create_maps();
399 404
      for (NodeIt it(*_gr); it != INVALID; ++it) {
400 405
	_pred->set(it, INVALID);
401 406
	_dist->set(it, value);
402 407
      }
403 408
      _process.clear();
404 409
      if (OperationTraits::less(value, OperationTraits::infinity())) {
405 410
	for (NodeIt it(*_gr); it != INVALID; ++it) {
406 411
	  _process.push_back(it);
407 412
	  _mask->set(it, true);
408 413
	}
409 414
      } else {
410 415
	for (NodeIt it(*_gr); it != INVALID; ++it) {
411 416
	  _mask->set(it, false);
412 417
	}
413 418
      }
414 419
    }
415 420
    
416 421
    /// \brief Adds a new source node.
417 422
    ///
418 423
    /// This function adds a new source node. The optional second parameter
419 424
    /// is the initial distance of the node.
420 425
    void addSource(Node source, Value dst = OperationTraits::zero()) {
421 426
      _dist->set(source, dst);
422 427
      if (!(*_mask)[source]) {
423 428
	_process.push_back(source);
424 429
	_mask->set(source, true);
425 430
      }
426 431
    }
427 432

	
428 433
    /// \brief Executes one round from the Bellman-Ford algorithm.
429 434
    ///
430 435
    /// If the algoritm calculated the distances in the previous round
431 436
    /// exactly for the paths of at most \c k arcs, then this function
432 437
    /// will calculate the distances exactly for the paths of at most
433 438
    /// <tt>k+1</tt> arcs. Performing \c k iterations using this function
434 439
    /// calculates the shortest path distances exactly for the paths
435 440
    /// consisting of at most \c k arcs.
436 441
    ///
437 442
    /// \warning The paths with limited arc number cannot be retrieved
438 443
    /// easily with \ref path() or \ref predArc() functions. If you also
439 444
    /// need the shortest paths and not only the distances, you should
440 445
    /// store the \ref predMap() "predecessor map" after each iteration
441 446
    /// and build the path manually.
442 447
    ///
443 448
    /// \return \c true when the algorithm have not found more shorter
444 449
    /// paths.
445 450
    ///
446 451
    /// \see ActiveIt
447 452
    bool processNextRound() {
448 453
      for (int i = 0; i < int(_process.size()); ++i) {
449 454
	_mask->set(_process[i], false);
450 455
      }
451 456
      std::vector<Node> nextProcess;
452 457
      std::vector<Value> values(_process.size());
453 458
      for (int i = 0; i < int(_process.size()); ++i) {
454 459
	values[i] = (*_dist)[_process[i]];
455 460
      }
456 461
      for (int i = 0; i < int(_process.size()); ++i) {
457 462
	for (OutArcIt it(*_gr, _process[i]); it != INVALID; ++it) {
458 463
	  Node target = _gr->target(it);
459 464
	  Value relaxed = OperationTraits::plus(values[i], (*_length)[it]);
460 465
	  if (OperationTraits::less(relaxed, (*_dist)[target])) {
461 466
	    _pred->set(target, it);
462 467
	    _dist->set(target, relaxed);
463 468
	    if (!(*_mask)[target]) {
464 469
	      _mask->set(target, true);
465 470
	      nextProcess.push_back(target);
466 471
	    }
467 472
	  }	  
468 473
	}
469 474
      }
470 475
      _process.swap(nextProcess);
471 476
      return _process.empty();
472 477
    }
473 478

	
474 479
    /// \brief Executes one weak round from the Bellman-Ford algorithm.
475 480
    ///
476 481
    /// If the algorithm calculated the distances in the previous round
477 482
    /// at least for the paths of at most \c k arcs, then this function
478 483
    /// will calculate the distances at least for the paths of at most
479 484
    /// <tt>k+1</tt> arcs.
480 485
    /// This function does not make it possible to calculate the shortest
481 486
    /// path distances exactly for paths consisting of at most \c k arcs,
482 487
    /// this is why it is called weak round.
483 488
    ///
484 489
    /// \return \c true when the algorithm have not found more shorter
485 490
    /// paths.
486 491
    ///
487 492
    /// \see ActiveIt
488 493
    bool processNextWeakRound() {
489 494
      for (int i = 0; i < int(_process.size()); ++i) {
490 495
	_mask->set(_process[i], false);
491 496
      }
492 497
      std::vector<Node> nextProcess;
493 498
      for (int i = 0; i < int(_process.size()); ++i) {
494 499
	for (OutArcIt it(*_gr, _process[i]); it != INVALID; ++it) {
495 500
	  Node target = _gr->target(it);
496 501
	  Value relaxed = 
497 502
	    OperationTraits::plus((*_dist)[_process[i]], (*_length)[it]);
498 503
	  if (OperationTraits::less(relaxed, (*_dist)[target])) {
499 504
	    _pred->set(target, it);
500 505
	    _dist->set(target, relaxed);
501 506
	    if (!(*_mask)[target]) {
502 507
	      _mask->set(target, true);
503 508
	      nextProcess.push_back(target);
504 509
	    }
505 510
	  }	  
506 511
	}
507 512
      }
508 513
      _process.swap(nextProcess);
509 514
      return _process.empty();
510 515
    }
511 516

	
512 517
    /// \brief Executes the algorithm.
513 518
    ///
514 519
    /// Executes the algorithm.
515 520
    ///
516 521
    /// This method runs the Bellman-Ford algorithm from the root node(s)
517 522
    /// in order to compute the shortest path to each node.
518 523
    ///
519 524
    /// The algorithm computes
520 525
    /// - the shortest path tree (forest),
521 526
    /// - the distance of each node from the root(s).
522 527
    ///
523 528
    /// \pre init() must be called and at least one root node should be
524 529
    /// added with addSource() before using this function.
525 530
    void start() {
526 531
      int num = countNodes(*_gr) - 1;
527 532
      for (int i = 0; i < num; ++i) {
528 533
	if (processNextWeakRound()) break;
529 534
      }
530 535
    }
531 536

	
532 537
    /// \brief Executes the algorithm and checks the negative cycles.
533 538
    ///
534 539
    /// Executes the algorithm and checks the negative cycles.
535 540
    ///
536 541
    /// This method runs the Bellman-Ford algorithm from the root node(s)
537 542
    /// in order to compute the shortest path to each node and also checks
538 543
    /// if the digraph contains cycles with negative total length.
539 544
    ///
540 545
    /// The algorithm computes 
541 546
    /// - the shortest path tree (forest),
542 547
    /// - the distance of each node from the root(s).
543 548
    /// 
544 549
    /// \return \c false if there is a negative cycle in the digraph.
545 550
    ///
546 551
    /// \pre init() must be called and at least one root node should be
547 552
    /// added with addSource() before using this function. 
548 553
    bool checkedStart() {
549 554
      int num = countNodes(*_gr);
550 555
      for (int i = 0; i < num; ++i) {
551 556
	if (processNextWeakRound()) return true;
552 557
      }
553 558
      return _process.empty();
554 559
    }
555 560

	
556 561
    /// \brief Executes the algorithm with arc number limit.
557 562
    ///
558 563
    /// Executes the algorithm with arc number limit.
559 564
    ///
560 565
    /// This method runs the Bellman-Ford algorithm from the root node(s)
561 566
    /// in order to compute the shortest path distance for each node
562 567
    /// using only the paths consisting of at most \c num arcs.
563 568
    ///
564 569
    /// The algorithm computes
565 570
    /// - the limited distance of each node from the root(s),
566 571
    /// - the predecessor arc for each node.
567 572
    ///
568 573
    /// \warning The paths with limited arc number cannot be retrieved
569 574
    /// easily with \ref path() or \ref predArc() functions. If you also
570 575
    /// need the shortest paths and not only the distances, you should
571 576
    /// store the \ref predMap() "predecessor map" after each iteration
572 577
    /// and build the path manually.
573 578
    ///
574 579
    /// \pre init() must be called and at least one root node should be
575 580
    /// added with addSource() before using this function. 
576 581
    void limitedStart(int num) {
577 582
      for (int i = 0; i < num; ++i) {
578 583
	if (processNextRound()) break;
579 584
      }
580 585
    }
581 586
    
582 587
    /// \brief Runs the algorithm from the given root node.
583 588
    ///    
584 589
    /// This method runs the Bellman-Ford algorithm from the given root
585 590
    /// node \c s in order to compute the shortest path to each node.
586 591
    ///
587 592
    /// The algorithm computes
588 593
    /// - the shortest path tree (forest),
589 594
    /// - the distance of each node from the root(s).
590 595
    ///
591 596
    /// \note bf.run(s) is just a shortcut of the following code.
592 597
    /// \code
593 598
    ///   bf.init();
594 599
    ///   bf.addSource(s);
595 600
    ///   bf.start();
596 601
    /// \endcode
597 602
    void run(Node s) {
598 603
      init();
599 604
      addSource(s);
600 605
      start();
601 606
    }
602 607
    
603 608
    /// \brief Runs the algorithm from the given root node with arc
604 609
    /// number limit.
605 610
    ///    
606 611
    /// This method runs the Bellman-Ford algorithm from the given root
607 612
    /// node \c s in order to compute the shortest path distance for each
608 613
    /// node using only the paths consisting of at most \c num arcs.
609 614
    ///
610 615
    /// The algorithm computes
611 616
    /// - the limited distance of each node from the root(s),
612 617
    /// - the predecessor arc for each node.
613 618
    ///
614 619
    /// \warning The paths with limited arc number cannot be retrieved
615 620
    /// easily with \ref path() or \ref predArc() functions. If you also
616 621
    /// need the shortest paths and not only the distances, you should
617 622
    /// store the \ref predMap() "predecessor map" after each iteration
618 623
    /// and build the path manually.
619 624
    ///
620 625
    /// \note bf.run(s, num) is just a shortcut of the following code.
621 626
    /// \code
622 627
    ///   bf.init();
623 628
    ///   bf.addSource(s);
624 629
    ///   bf.limitedStart(num);
625 630
    /// \endcode
626 631
    void run(Node s, int num) {
627 632
      init();
628 633
      addSource(s);
629 634
      limitedStart(num);
630 635
    }
631 636
    
632 637
    ///@}
633 638

	
634 639
    /// \brief LEMON iterator for getting the active nodes.
635 640
    ///
636 641
    /// This class provides a common style LEMON iterator that traverses
637 642
    /// the active nodes of the Bellman-Ford algorithm after the last
638 643
    /// phase. These nodes should be checked in the next phase to
639 644
    /// find augmenting arcs outgoing from them.
640 645
    class ActiveIt {
641 646
    public:
642 647

	
643 648
      /// \brief Constructor.
644 649
      ///
645 650
      /// Constructor for getting the active nodes of the given BellmanFord
646 651
      /// instance. 
647 652
      ActiveIt(const BellmanFord& algorithm) : _algorithm(&algorithm)
648 653
      {
649 654
        _index = _algorithm->_process.size() - 1;
650 655
      }
651 656

	
652 657
      /// \brief Invalid constructor.
653 658
      ///
654 659
      /// Invalid constructor.
655 660
      ActiveIt(Invalid) : _algorithm(0), _index(-1) {}
656 661

	
657 662
      /// \brief Conversion to \c Node.
658 663
      ///
659 664
      /// Conversion to \c Node.
660 665
      operator Node() const { 
661 666
        return _index >= 0 ? _algorithm->_process[_index] : INVALID;
662 667
      }
663 668

	
664 669
      /// \brief Increment operator.
665 670
      ///
666 671
      /// Increment operator.
667 672
      ActiveIt& operator++() {
668 673
        --_index;
669 674
        return *this; 
670 675
      }
671 676

	
672 677
      bool operator==(const ActiveIt& it) const { 
673 678
        return static_cast<Node>(*this) == static_cast<Node>(it); 
674 679
      }
675 680
      bool operator!=(const ActiveIt& it) const { 
676 681
        return static_cast<Node>(*this) != static_cast<Node>(it); 
677 682
      }
678 683
      bool operator<(const ActiveIt& it) const { 
679 684
        return static_cast<Node>(*this) < static_cast<Node>(it); 
680 685
      }
681 686
      
682 687
    private:
683 688
      const BellmanFord* _algorithm;
684 689
      int _index;
685 690
    };
686 691
    
687 692
    /// \name Query Functions
688 693
    /// The result of the Bellman-Ford algorithm can be obtained using these
689 694
    /// functions.\n
690 695
    /// Either \ref run() or \ref init() should be called before using them.
691 696
    
692 697
    ///@{
693 698

	
694 699
    /// \brief The shortest path to the given node.
695 700
    ///    
696 701
    /// Gives back the shortest path to the given node from the root(s).
697 702
    ///
698 703
    /// \warning \c t should be reached from the root(s).
699 704
    ///
700 705
    /// \pre Either \ref run() or \ref init() must be called before
701 706
    /// using this function.
702 707
    Path path(Node t) const
703 708
    {
704 709
      return Path(*_gr, *_pred, t);
705 710
    }
706 711
	  
707 712
    /// \brief The distance of the given node from the root(s).
708 713
    ///
709 714
    /// Returns the distance of the given node from the root(s).
710 715
    ///
711 716
    /// \warning If node \c v is not reached from the root(s), then
712 717
    /// the return value of this function is undefined.
713 718
    ///
714 719
    /// \pre Either \ref run() or \ref init() must be called before
715 720
    /// using this function.
716 721
    Value dist(Node v) const { return (*_dist)[v]; }
717 722

	
718 723
    /// \brief Returns the 'previous arc' of the shortest path tree for
719 724
    /// the given node.
720 725
    ///
721 726
    /// This function returns the 'previous arc' of the shortest path
722 727
    /// tree for node \c v, i.e. it returns the last arc of a
723 728
    /// shortest path from a root to \c v. It is \c INVALID if \c v
724 729
    /// is not reached from the root(s) or if \c v is a root.
725 730
    ///
726 731
    /// The shortest path tree used here is equal to the shortest path
727 732
    /// tree used in \ref predNode() and \ref predMap().
728 733
    ///
729 734
    /// \pre Either \ref run() or \ref init() must be called before
730 735
    /// using this function.
731 736
    Arc predArc(Node v) const { return (*_pred)[v]; }
732 737

	
733 738
    /// \brief Returns the 'previous node' of the shortest path tree for
734 739
    /// the given node.
735 740
    ///
736 741
    /// This function returns the 'previous node' of the shortest path
737 742
    /// tree for node \c v, i.e. it returns the last but one node of
738 743
    /// a shortest path from a root to \c v. It is \c INVALID if \c v
739 744
    /// is not reached from the root(s) or if \c v is a root.
740 745
    ///
741 746
    /// The shortest path tree used here is equal to the shortest path
742 747
    /// tree used in \ref predArc() and \ref predMap().
743 748
    ///
744 749
    /// \pre Either \ref run() or \ref init() must be called before
745 750
    /// using this function.
746 751
    Node predNode(Node v) const { 
747 752
      return (*_pred)[v] == INVALID ? INVALID : _gr->source((*_pred)[v]); 
748 753
    }
749 754
    
750 755
    /// \brief Returns a const reference to the node map that stores the
751 756
    /// distances of the nodes.
752 757
    ///
753 758
    /// Returns a const reference to the node map that stores the distances
754 759
    /// of the nodes calculated by the algorithm.
755 760
    ///
756 761
    /// \pre Either \ref run() or \ref init() must be called before
757 762
    /// using this function.
758 763
    const DistMap &distMap() const { return *_dist;}
759 764
 
760 765
    /// \brief Returns a const reference to the node map that stores the
761 766
    /// predecessor arcs.
762 767
    ///
763 768
    /// Returns a const reference to the node map that stores the predecessor
764 769
    /// arcs, which form the shortest path tree (forest).
765 770
    ///
766 771
    /// \pre Either \ref run() or \ref init() must be called before
767 772
    /// using this function.
768 773
    const PredMap &predMap() const { return *_pred; }
769 774
 
770 775
    /// \brief Checks if a node is reached from the root(s).
771 776
    ///
772 777
    /// Returns \c true if \c v is reached from the root(s).
773 778
    ///
774 779
    /// \pre Either \ref run() or \ref init() must be called before
775 780
    /// using this function.
776 781
    bool reached(Node v) const {
777 782
      return (*_dist)[v] != OperationTraits::infinity();
778 783
    }
779 784

	
780 785
    /// \brief Gives back a negative cycle.
781 786
    ///    
782 787
    /// This function gives back a directed cycle with negative total
783 788
    /// length if the algorithm has already found one.
784 789
    /// Otherwise it gives back an empty path.
785 790
    lemon::Path<Digraph> negativeCycle() const {
786 791
      typename Digraph::template NodeMap<int> state(*_gr, -1);
787 792
      lemon::Path<Digraph> cycle;
788 793
      for (int i = 0; i < int(_process.size()); ++i) {
789 794
        if (state[_process[i]] != -1) continue;
790 795
        for (Node v = _process[i]; (*_pred)[v] != INVALID;
791 796
             v = _gr->source((*_pred)[v])) {
792 797
          if (state[v] == i) {
793 798
            cycle.addFront((*_pred)[v]);
794 799
            for (Node u = _gr->source((*_pred)[v]); u != v;
795 800
                 u = _gr->source((*_pred)[u])) {
796 801
              cycle.addFront((*_pred)[u]);
797 802
            }
798 803
            return cycle;
799 804
          }
800 805
          else if (state[v] >= 0) {
801 806
            break;
802 807
          }
803 808
          state[v] = i;
804 809
        }
805 810
      }
806 811
      return cycle;
807 812
    }
808 813
    
809 814
    ///@}
810 815
  };
811 816
 
812 817
  /// \brief Default traits class of bellmanFord() function.
813 818
  ///
814 819
  /// Default traits class of bellmanFord() function.
815 820
  /// \tparam GR The type of the digraph.
816 821
  /// \tparam LEN The type of the length map.
817 822
  template <typename GR, typename LEN>
818 823
  struct BellmanFordWizardDefaultTraits {
819 824
    /// The type of the digraph the algorithm runs on. 
820 825
    typedef GR Digraph;
821 826

	
822 827
    /// \brief The type of the map that stores the arc lengths.
823 828
    ///
824 829
    /// The type of the map that stores the arc lengths.
825 830
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
826 831
    typedef LEN LengthMap;
827 832

	
828 833
    /// The type of the arc lengths.
829 834
    typedef typename LEN::Value Value;
830 835

	
831 836
    /// \brief Operation traits for Bellman-Ford algorithm.
832 837
    ///
833 838
    /// It defines the used operations and the infinity value for the
834 839
    /// given \c Value type.
835 840
    /// \see BellmanFordDefaultOperationTraits
836 841
    typedef BellmanFordDefaultOperationTraits<Value> OperationTraits;
837 842

	
838 843
    /// \brief The type of the map that stores the last
839 844
    /// arcs of the shortest paths.
840 845
    /// 
841 846
    /// The type of the map that stores the last arcs of the shortest paths.
842 847
    /// It must conform to the \ref concepts::WriteMap "WriteMap" concept.
843 848
    typedef typename GR::template NodeMap<typename GR::Arc> PredMap;
844 849

	
845 850
    /// \brief Instantiates a \c PredMap.
846 851
    /// 
847 852
    /// This function instantiates a \ref PredMap.
848 853
    /// \param g is the digraph to which we would like to define the
849 854
    /// \ref PredMap.
850 855
    static PredMap *createPredMap(const GR &g) {
851 856
      return new PredMap(g);
852 857
    }
853 858

	
854 859
    /// \brief The type of the map that stores the distances of the nodes.
855 860
    ///
856 861
    /// The type of the map that stores the distances of the nodes.
857 862
    /// It must conform to the \ref concepts::WriteMap "WriteMap" concept.
858 863
    typedef typename GR::template NodeMap<Value> DistMap;
859 864

	
860 865
    /// \brief Instantiates a \c DistMap.
861 866
    ///
862 867
    /// This function instantiates a \ref DistMap. 
863 868
    /// \param g is the digraph to which we would like to define the
864 869
    /// \ref DistMap.
865 870
    static DistMap *createDistMap(const GR &g) {
866 871
      return new DistMap(g);
867 872
    }
868 873

	
869 874
    ///The type of the shortest paths.
870 875

	
871 876
    ///The type of the shortest paths.
872 877
    ///It must meet the \ref concepts::Path "Path" concept.
873 878
    typedef lemon::Path<Digraph> Path;
874 879
  };
875 880
  
876 881
  /// \brief Default traits class used by BellmanFordWizard.
877 882
  ///
878 883
  /// Default traits class used by BellmanFordWizard.
879 884
  /// \tparam GR The type of the digraph.
880 885
  /// \tparam LEN The type of the length map.
881 886
  template <typename GR, typename LEN>
882 887
  class BellmanFordWizardBase 
883 888
    : public BellmanFordWizardDefaultTraits<GR, LEN> {
884 889

	
885 890
    typedef BellmanFordWizardDefaultTraits<GR, LEN> Base;
886 891
  protected:
887 892
    // Type of the nodes in the digraph.
888 893
    typedef typename Base::Digraph::Node Node;
889 894

	
890 895
    // Pointer to the underlying digraph.
891 896
    void *_graph;
892 897
    // Pointer to the length map
893 898
    void *_length;
894 899
    // Pointer to the map of predecessors arcs.
895 900
    void *_pred;
896 901
    // Pointer to the map of distances.
897 902
    void *_dist;
898 903
    //Pointer to the shortest path to the target node.
899 904
    void *_path;
900 905
    //Pointer to the distance of the target node.
901 906
    void *_di;
902 907

	
903 908
    public:
904 909
    /// Constructor.
905 910
    
906 911
    /// This constructor does not require parameters, it initiates
907 912
    /// all of the attributes to default values \c 0.
908 913
    BellmanFordWizardBase() :
909 914
      _graph(0), _length(0), _pred(0), _dist(0), _path(0), _di(0) {}
910 915

	
911 916
    /// Constructor.
912 917
    
913 918
    /// This constructor requires two parameters,
914 919
    /// others are initiated to \c 0.
915 920
    /// \param gr The digraph the algorithm runs on.
916 921
    /// \param len The length map.
917 922
    BellmanFordWizardBase(const GR& gr, 
918 923
			  const LEN& len) :
919 924
      _graph(reinterpret_cast<void*>(const_cast<GR*>(&gr))), 
920 925
      _length(reinterpret_cast<void*>(const_cast<LEN*>(&len))), 
921 926
      _pred(0), _dist(0), _path(0), _di(0) {}
922 927

	
923 928
  };
924 929
  
925 930
  /// \brief Auxiliary class for the function-type interface of the
926 931
  /// \ref BellmanFord "Bellman-Ford" algorithm.
927 932
  ///
928 933
  /// This auxiliary class is created to implement the
929 934
  /// \ref bellmanFord() "function-type interface" of the
930 935
  /// \ref BellmanFord "Bellman-Ford" algorithm.
931 936
  /// It does not have own \ref run() method, it uses the
932 937
  /// functions and features of the plain \ref BellmanFord.
933 938
  ///
934 939
  /// This class should only be used through the \ref bellmanFord()
935 940
  /// function, which makes it easier to use the algorithm.
941
  ///
942
  /// \tparam TR The traits class that defines various types used by the
943
  /// algorithm.
936 944
  template<class TR>
937 945
  class BellmanFordWizard : public TR {
938 946
    typedef TR Base;
939 947

	
940 948
    typedef typename TR::Digraph Digraph;
941 949

	
942 950
    typedef typename Digraph::Node Node;
943 951
    typedef typename Digraph::NodeIt NodeIt;
944 952
    typedef typename Digraph::Arc Arc;
945 953
    typedef typename Digraph::OutArcIt ArcIt;
946 954
    
947 955
    typedef typename TR::LengthMap LengthMap;
948 956
    typedef typename LengthMap::Value Value;
949 957
    typedef typename TR::PredMap PredMap;
950 958
    typedef typename TR::DistMap DistMap;
951 959
    typedef typename TR::Path Path;
952 960

	
953 961
  public:
954 962
    /// Constructor.
955 963
    BellmanFordWizard() : TR() {}
956 964

	
957 965
    /// \brief Constructor that requires parameters.
958 966
    ///
959 967
    /// Constructor that requires parameters.
960 968
    /// These parameters will be the default values for the traits class.
961 969
    /// \param gr The digraph the algorithm runs on.
962 970
    /// \param len The length map.
963 971
    BellmanFordWizard(const Digraph& gr, const LengthMap& len) 
964 972
      : TR(gr, len) {}
965 973

	
966 974
    /// \brief Copy constructor
967 975
    BellmanFordWizard(const TR &b) : TR(b) {}
968 976

	
969 977
    ~BellmanFordWizard() {}
970 978

	
971 979
    /// \brief Runs the Bellman-Ford algorithm from the given source node.
972 980
    ///    
973 981
    /// This method runs the Bellman-Ford algorithm from the given source
974 982
    /// node in order to compute the shortest path to each node.
975 983
    void run(Node s) {
976 984
      BellmanFord<Digraph,LengthMap,TR> 
977 985
	bf(*reinterpret_cast<const Digraph*>(Base::_graph), 
978 986
           *reinterpret_cast<const LengthMap*>(Base::_length));
979 987
      if (Base::_pred) bf.predMap(*reinterpret_cast<PredMap*>(Base::_pred));
980 988
      if (Base::_dist) bf.distMap(*reinterpret_cast<DistMap*>(Base::_dist));
981 989
      bf.run(s);
982 990
    }
983 991

	
984 992
    /// \brief Runs the Bellman-Ford algorithm to find the shortest path
985 993
    /// between \c s and \c t.
986 994
    ///
987 995
    /// This method runs the Bellman-Ford algorithm from node \c s
988 996
    /// in order to compute the shortest path to node \c t.
989 997
    /// Actually, it computes the shortest path to each node, but using
990 998
    /// this function you can retrieve the distance and the shortest path
991 999
    /// for a single target node easier.
992 1000
    ///
993 1001
    /// \return \c true if \c t is reachable form \c s.
994 1002
    bool run(Node s, Node t) {
995 1003
      BellmanFord<Digraph,LengthMap,TR>
996 1004
        bf(*reinterpret_cast<const Digraph*>(Base::_graph),
997 1005
           *reinterpret_cast<const LengthMap*>(Base::_length));
998 1006
      if (Base::_pred) bf.predMap(*reinterpret_cast<PredMap*>(Base::_pred));
999 1007
      if (Base::_dist) bf.distMap(*reinterpret_cast<DistMap*>(Base::_dist));
1000 1008
      bf.run(s);
1001 1009
      if (Base::_path) *reinterpret_cast<Path*>(Base::_path) = bf.path(t);
1002 1010
      if (Base::_di) *reinterpret_cast<Value*>(Base::_di) = bf.dist(t);
1003 1011
      return bf.reached(t);
1004 1012
    }
1005 1013

	
1006 1014
    template<class T>
1007 1015
    struct SetPredMapBase : public Base {
1008 1016
      typedef T PredMap;
1009 1017
      static PredMap *createPredMap(const Digraph &) { return 0; };
1010 1018
      SetPredMapBase(const TR &b) : TR(b) {}
1011 1019
    };
1012 1020
    
1013 1021
    /// \brief \ref named-templ-param "Named parameter" for setting
1014 1022
    /// the predecessor map.
1015 1023
    ///
1016 1024
    /// \ref named-templ-param "Named parameter" for setting
1017 1025
    /// the map that stores the predecessor arcs of the nodes.
1018 1026
    template<class T>
1019 1027
    BellmanFordWizard<SetPredMapBase<T> > predMap(const T &t) {
1020 1028
      Base::_pred=reinterpret_cast<void*>(const_cast<T*>(&t));
1021 1029
      return BellmanFordWizard<SetPredMapBase<T> >(*this);
1022 1030
    }
1023 1031
    
1024 1032
    template<class T>
1025 1033
    struct SetDistMapBase : public Base {
1026 1034
      typedef T DistMap;
1027 1035
      static DistMap *createDistMap(const Digraph &) { return 0; };
1028 1036
      SetDistMapBase(const TR &b) : TR(b) {}
1029 1037
    };
1030 1038
    
1031 1039
    /// \brief \ref named-templ-param "Named parameter" for setting
1032 1040
    /// the distance map.
1033 1041
    ///
1034 1042
    /// \ref named-templ-param "Named parameter" for setting
1035 1043
    /// the map that stores the distances of the nodes calculated
1036 1044
    /// by the algorithm.
1037 1045
    template<class T>
1038 1046
    BellmanFordWizard<SetDistMapBase<T> > distMap(const T &t) {
1039 1047
      Base::_dist=reinterpret_cast<void*>(const_cast<T*>(&t));
1040 1048
      return BellmanFordWizard<SetDistMapBase<T> >(*this);
1041 1049
    }
1042 1050

	
1043 1051
    template<class T>
1044 1052
    struct SetPathBase : public Base {
1045 1053
      typedef T Path;
1046 1054
      SetPathBase(const TR &b) : TR(b) {}
1047 1055
    };
1048 1056

	
1049 1057
    /// \brief \ref named-func-param "Named parameter" for getting
1050 1058
    /// the shortest path to the target node.
1051 1059
    ///
1052 1060
    /// \ref named-func-param "Named parameter" for getting
1053 1061
    /// the shortest path to the target node.
1054 1062
    template<class T>
1055 1063
    BellmanFordWizard<SetPathBase<T> > path(const T &t)
1056 1064
    {
1057 1065
      Base::_path=reinterpret_cast<void*>(const_cast<T*>(&t));
1058 1066
      return BellmanFordWizard<SetPathBase<T> >(*this);
1059 1067
    }
1060 1068

	
1061 1069
    /// \brief \ref named-func-param "Named parameter" for getting
1062 1070
    /// the distance of the target node.
1063 1071
    ///
1064 1072
    /// \ref named-func-param "Named parameter" for getting
1065 1073
    /// the distance of the target node.
1066 1074
    BellmanFordWizard dist(const Value &d)
1067 1075
    {
1068 1076
      Base::_di=reinterpret_cast<void*>(const_cast<Value*>(&d));
1069 1077
      return *this;
1070 1078
    }
1071 1079
    
1072 1080
  };
1073 1081
  
1074 1082
  /// \brief Function type interface for the \ref BellmanFord "Bellman-Ford"
1075 1083
  /// algorithm.
1076 1084
  ///
1077 1085
  /// \ingroup shortest_path
1078 1086
  /// Function type interface for the \ref BellmanFord "Bellman-Ford"
1079 1087
  /// algorithm.
1080 1088
  ///
1081 1089
  /// This function also has several \ref named-templ-func-param 
1082 1090
  /// "named parameters", they are declared as the members of class 
1083 1091
  /// \ref BellmanFordWizard.
1084 1092
  /// The following examples show how to use these parameters.
1085 1093
  /// \code
1086 1094
  ///   // Compute shortest path from node s to each node
1087 1095
  ///   bellmanFord(g,length).predMap(preds).distMap(dists).run(s);
1088 1096
  ///
1089 1097
  ///   // Compute shortest path from s to t
1090 1098
  ///   bool reached = bellmanFord(g,length).path(p).dist(d).run(s,t);
1091 1099
  /// \endcode
1092 1100
  /// \warning Don't forget to put the \ref BellmanFordWizard::run() "run()"
1093 1101
  /// to the end of the parameter list.
1094 1102
  /// \sa BellmanFordWizard
1095 1103
  /// \sa BellmanFord
1096 1104
  template<typename GR, typename LEN>
1097 1105
  BellmanFordWizard<BellmanFordWizardBase<GR,LEN> >
1098 1106
  bellmanFord(const GR& digraph,
1099 1107
	      const LEN& length)
1100 1108
  {
1101 1109
    return BellmanFordWizard<BellmanFordWizardBase<GR,LEN> >(digraph, length);
1102 1110
  }
1103 1111

	
1104 1112
} //END OF NAMESPACE LEMON
1105 1113

	
1106 1114
#endif
1107 1115

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

	
19 19
#ifndef LEMON_BFS_H
20 20
#define LEMON_BFS_H
21 21

	
22 22
///\ingroup search
23 23
///\file
24 24
///\brief BFS algorithm.
25 25

	
26 26
#include <lemon/list_graph.h>
27 27
#include <lemon/bits/path_dump.h>
28 28
#include <lemon/core.h>
29 29
#include <lemon/error.h>
30 30
#include <lemon/maps.h>
31 31
#include <lemon/path.h>
32 32

	
33 33
namespace lemon {
34 34

	
35 35
  ///Default traits class of Bfs class.
36 36

	
37 37
  ///Default traits class of Bfs class.
38 38
  ///\tparam GR Digraph type.
39 39
  template<class GR>
40 40
  struct BfsDefaultTraits
41 41
  {
42 42
    ///The type of the digraph the algorithm runs on.
43 43
    typedef GR Digraph;
44 44

	
45 45
    ///\brief The type of the map that stores the predecessor
46 46
    ///arcs of the shortest paths.
47 47
    ///
48 48
    ///The type of the map that stores the predecessor
49 49
    ///arcs of the shortest paths.
50 50
    ///It must conform to the \ref concepts::WriteMap "WriteMap" concept.
51 51
    typedef typename Digraph::template NodeMap<typename Digraph::Arc> PredMap;
52 52
    ///Instantiates a \c PredMap.
53 53

	
54 54
    ///This function instantiates a \ref PredMap.
55 55
    ///\param g is the digraph, to which we would like to define the
56 56
    ///\ref PredMap.
57 57
    static PredMap *createPredMap(const Digraph &g)
58 58
    {
59 59
      return new PredMap(g);
60 60
    }
61 61

	
62 62
    ///The type of the map that indicates which nodes are processed.
63 63

	
64 64
    ///The type of the map that indicates which nodes are processed.
65 65
    ///It must conform to the \ref concepts::WriteMap "WriteMap" concept.
66 66
    ///By default, it is a NullMap.
67 67
    typedef NullMap<typename Digraph::Node,bool> ProcessedMap;
68 68
    ///Instantiates a \c ProcessedMap.
69 69

	
70 70
    ///This function instantiates a \ref ProcessedMap.
71 71
    ///\param g is the digraph, to which
72 72
    ///we would like to define the \ref ProcessedMap
73 73
#ifdef DOXYGEN
74 74
    static ProcessedMap *createProcessedMap(const Digraph &g)
75 75
#else
76 76
    static ProcessedMap *createProcessedMap(const Digraph &)
77 77
#endif
78 78
    {
79 79
      return new ProcessedMap();
80 80
    }
81 81

	
82 82
    ///The type of the map that indicates which nodes are reached.
83 83

	
84 84
    ///The type of the map that indicates which nodes are reached.
85 85
    ///It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
86 86
    typedef typename Digraph::template NodeMap<bool> ReachedMap;
87 87
    ///Instantiates a \c ReachedMap.
88 88

	
89 89
    ///This function instantiates a \ref ReachedMap.
90 90
    ///\param g is the digraph, to which
91 91
    ///we would like to define the \ref ReachedMap.
92 92
    static ReachedMap *createReachedMap(const Digraph &g)
93 93
    {
94 94
      return new ReachedMap(g);
95 95
    }
96 96

	
97 97
    ///The type of the map that stores the distances of the nodes.
98 98

	
99 99
    ///The type of the map that stores the distances of the nodes.
100 100
    ///It must conform to the \ref concepts::WriteMap "WriteMap" concept.
101 101
    typedef typename Digraph::template NodeMap<int> DistMap;
102 102
    ///Instantiates a \c DistMap.
103 103

	
104 104
    ///This function instantiates a \ref DistMap.
105 105
    ///\param g is the digraph, to which we would like to define the
106 106
    ///\ref DistMap.
107 107
    static DistMap *createDistMap(const Digraph &g)
108 108
    {
109 109
      return new DistMap(g);
110 110
    }
111 111
  };
112 112

	
113 113
  ///%BFS algorithm class.
114 114

	
115 115
  ///\ingroup search
116 116
  ///This class provides an efficient implementation of the %BFS algorithm.
117 117
  ///
118 118
  ///There is also a \ref bfs() "function-type interface" for the BFS
119 119
  ///algorithm, which is convenient in the simplier cases and it can be
120 120
  ///used easier.
121 121
  ///
122 122
  ///\tparam GR The type of the digraph the algorithm runs on.
123 123
  ///The default type is \ref ListDigraph.
124
  ///\tparam TR The traits class that defines various types used by the
125
  ///algorithm. By default, it is \ref BfsDefaultTraits
126
  ///"BfsDefaultTraits<GR>".
127
  ///In most cases, this parameter should not be set directly,
128
  ///consider to use the named template parameters instead.
124 129
#ifdef DOXYGEN
125 130
  template <typename GR,
126 131
            typename TR>
127 132
#else
128 133
  template <typename GR=ListDigraph,
129 134
            typename TR=BfsDefaultTraits<GR> >
130 135
#endif
131 136
  class Bfs {
132 137
  public:
133 138

	
134 139
    ///The type of the digraph the algorithm runs on.
135 140
    typedef typename TR::Digraph Digraph;
136 141

	
137 142
    ///\brief The type of the map that stores the predecessor arcs of the
138 143
    ///shortest paths.
139 144
    typedef typename TR::PredMap PredMap;
140 145
    ///The type of the map that stores the distances of the nodes.
141 146
    typedef typename TR::DistMap DistMap;
142 147
    ///The type of the map that indicates which nodes are reached.
143 148
    typedef typename TR::ReachedMap ReachedMap;
144 149
    ///The type of the map that indicates which nodes are processed.
145 150
    typedef typename TR::ProcessedMap ProcessedMap;
146 151
    ///The type of the paths.
147 152
    typedef PredMapPath<Digraph, PredMap> Path;
148 153

	
149 154
    ///The \ref BfsDefaultTraits "traits class" of the algorithm.
150 155
    typedef TR Traits;
151 156

	
152 157
  private:
153 158

	
154 159
    typedef typename Digraph::Node Node;
155 160
    typedef typename Digraph::NodeIt NodeIt;
156 161
    typedef typename Digraph::Arc Arc;
157 162
    typedef typename Digraph::OutArcIt OutArcIt;
158 163

	
159 164
    //Pointer to the underlying digraph.
160 165
    const Digraph *G;
161 166
    //Pointer to the map of predecessor arcs.
162 167
    PredMap *_pred;
163 168
    //Indicates if _pred is locally allocated (true) or not.
164 169
    bool local_pred;
165 170
    //Pointer to the map of distances.
166 171
    DistMap *_dist;
167 172
    //Indicates if _dist is locally allocated (true) or not.
168 173
    bool local_dist;
169 174
    //Pointer to the map of reached status of the nodes.
170 175
    ReachedMap *_reached;
171 176
    //Indicates if _reached is locally allocated (true) or not.
172 177
    bool local_reached;
173 178
    //Pointer to the map of processed status of the nodes.
174 179
    ProcessedMap *_processed;
175 180
    //Indicates if _processed is locally allocated (true) or not.
176 181
    bool local_processed;
177 182

	
178 183
    std::vector<typename Digraph::Node> _queue;
179 184
    int _queue_head,_queue_tail,_queue_next_dist;
180 185
    int _curr_dist;
181 186

	
182 187
    //Creates the maps if necessary.
183 188
    void create_maps()
184 189
    {
185 190
      if(!_pred) {
186 191
        local_pred = true;
187 192
        _pred = Traits::createPredMap(*G);
188 193
      }
189 194
      if(!_dist) {
190 195
        local_dist = true;
191 196
        _dist = Traits::createDistMap(*G);
192 197
      }
193 198
      if(!_reached) {
194 199
        local_reached = true;
195 200
        _reached = Traits::createReachedMap(*G);
196 201
      }
197 202
      if(!_processed) {
198 203
        local_processed = true;
199 204
        _processed = Traits::createProcessedMap(*G);
200 205
      }
201 206
    }
202 207

	
203 208
  protected:
204 209

	
205 210
    Bfs() {}
206 211

	
207 212
  public:
208 213

	
209 214
    typedef Bfs Create;
210 215

	
211 216
    ///\name Named Template Parameters
212 217

	
213 218
    ///@{
214 219

	
215 220
    template <class T>
216 221
    struct SetPredMapTraits : public Traits {
217 222
      typedef T PredMap;
218 223
      static PredMap *createPredMap(const Digraph &)
219 224
      {
220 225
        LEMON_ASSERT(false, "PredMap is not initialized");
221 226
        return 0; // ignore warnings
222 227
      }
223 228
    };
224 229
    ///\brief \ref named-templ-param "Named parameter" for setting
225 230
    ///\c PredMap type.
226 231
    ///
227 232
    ///\ref named-templ-param "Named parameter" for setting
228 233
    ///\c PredMap type.
229 234
    ///It must conform to the \ref concepts::WriteMap "WriteMap" concept.
230 235
    template <class T>
231 236
    struct SetPredMap : public Bfs< Digraph, SetPredMapTraits<T> > {
232 237
      typedef Bfs< Digraph, SetPredMapTraits<T> > Create;
233 238
    };
234 239

	
235 240
    template <class T>
236 241
    struct SetDistMapTraits : public Traits {
237 242
      typedef T DistMap;
238 243
      static DistMap *createDistMap(const Digraph &)
239 244
      {
240 245
        LEMON_ASSERT(false, "DistMap is not initialized");
241 246
        return 0; // ignore warnings
242 247
      }
243 248
    };
244 249
    ///\brief \ref named-templ-param "Named parameter" for setting
245 250
    ///\c DistMap type.
246 251
    ///
247 252
    ///\ref named-templ-param "Named parameter" for setting
248 253
    ///\c DistMap type.
249 254
    ///It must conform to the \ref concepts::WriteMap "WriteMap" concept.
250 255
    template <class T>
251 256
    struct SetDistMap : public Bfs< Digraph, SetDistMapTraits<T> > {
252 257
      typedef Bfs< Digraph, SetDistMapTraits<T> > Create;
253 258
    };
254 259

	
255 260
    template <class T>
256 261
    struct SetReachedMapTraits : public Traits {
257 262
      typedef T ReachedMap;
258 263
      static ReachedMap *createReachedMap(const Digraph &)
259 264
      {
260 265
        LEMON_ASSERT(false, "ReachedMap is not initialized");
261 266
        return 0; // ignore warnings
262 267
      }
263 268
    };
264 269
    ///\brief \ref named-templ-param "Named parameter" for setting
265 270
    ///\c ReachedMap type.
266 271
    ///
267 272
    ///\ref named-templ-param "Named parameter" for setting
268 273
    ///\c ReachedMap type.
269 274
    ///It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
270 275
    template <class T>
271 276
    struct SetReachedMap : public Bfs< Digraph, SetReachedMapTraits<T> > {
272 277
      typedef Bfs< Digraph, SetReachedMapTraits<T> > Create;
273 278
    };
274 279

	
275 280
    template <class T>
276 281
    struct SetProcessedMapTraits : public Traits {
277 282
      typedef T ProcessedMap;
278 283
      static ProcessedMap *createProcessedMap(const Digraph &)
279 284
      {
280 285
        LEMON_ASSERT(false, "ProcessedMap is not initialized");
281 286
        return 0; // ignore warnings
282 287
      }
283 288
    };
284 289
    ///\brief \ref named-templ-param "Named parameter" for setting
285 290
    ///\c ProcessedMap type.
286 291
    ///
287 292
    ///\ref named-templ-param "Named parameter" for setting
288 293
    ///\c ProcessedMap type.
289 294
    ///It must conform to the \ref concepts::WriteMap "WriteMap" concept.
290 295
    template <class T>
291 296
    struct SetProcessedMap : public Bfs< Digraph, SetProcessedMapTraits<T> > {
292 297
      typedef Bfs< Digraph, SetProcessedMapTraits<T> > Create;
293 298
    };
294 299

	
295 300
    struct SetStandardProcessedMapTraits : public Traits {
296 301
      typedef typename Digraph::template NodeMap<bool> ProcessedMap;
297 302
      static ProcessedMap *createProcessedMap(const Digraph &g)
298 303
      {
299 304
        return new ProcessedMap(g);
300 305
        return 0; // ignore warnings
301 306
      }
302 307
    };
303 308
    ///\brief \ref named-templ-param "Named parameter" for setting
304 309
    ///\c ProcessedMap type to be <tt>Digraph::NodeMap<bool></tt>.
305 310
    ///
306 311
    ///\ref named-templ-param "Named parameter" for setting
307 312
    ///\c ProcessedMap type to be <tt>Digraph::NodeMap<bool></tt>.
308 313
    ///If you don't set it explicitly, it will be automatically allocated.
309 314
    struct SetStandardProcessedMap :
310 315
      public Bfs< Digraph, SetStandardProcessedMapTraits > {
311 316
      typedef Bfs< Digraph, SetStandardProcessedMapTraits > Create;
312 317
    };
313 318

	
314 319
    ///@}
315 320

	
316 321
  public:
317 322

	
318 323
    ///Constructor.
319 324

	
320 325
    ///Constructor.
321 326
    ///\param g The digraph the algorithm runs on.
322 327
    Bfs(const Digraph &g) :
323 328
      G(&g),
324 329
      _pred(NULL), local_pred(false),
325 330
      _dist(NULL), local_dist(false),
326 331
      _reached(NULL), local_reached(false),
327 332
      _processed(NULL), local_processed(false)
328 333
    { }
329 334

	
330 335
    ///Destructor.
331 336
    ~Bfs()
332 337
    {
333 338
      if(local_pred) delete _pred;
334 339
      if(local_dist) delete _dist;
335 340
      if(local_reached) delete _reached;
336 341
      if(local_processed) delete _processed;
337 342
    }
338 343

	
339 344
    ///Sets the map that stores the predecessor arcs.
340 345

	
341 346
    ///Sets the map that stores the predecessor arcs.
342 347
    ///If you don't use this function before calling \ref run(Node) "run()"
343 348
    ///or \ref init(), an instance will be allocated automatically.
344 349
    ///The destructor deallocates this automatically allocated map,
345 350
    ///of course.
346 351
    ///\return <tt> (*this) </tt>
347 352
    Bfs &predMap(PredMap &m)
348 353
    {
349 354
      if(local_pred) {
350 355
        delete _pred;
351 356
        local_pred=false;
352 357
      }
353 358
      _pred = &m;
354 359
      return *this;
355 360
    }
356 361

	
357 362
    ///Sets the map that indicates which nodes are reached.
358 363

	
359 364
    ///Sets the map that indicates which nodes are reached.
360 365
    ///If you don't use this function before calling \ref run(Node) "run()"
361 366
    ///or \ref init(), an instance will be allocated automatically.
362 367
    ///The destructor deallocates this automatically allocated map,
363 368
    ///of course.
364 369
    ///\return <tt> (*this) </tt>
365 370
    Bfs &reachedMap(ReachedMap &m)
366 371
    {
367 372
      if(local_reached) {
368 373
        delete _reached;
369 374
        local_reached=false;
370 375
      }
371 376
      _reached = &m;
372 377
      return *this;
373 378
    }
374 379

	
375 380
    ///Sets the map that indicates which nodes are processed.
376 381

	
377 382
    ///Sets the map that indicates which nodes are processed.
378 383
    ///If you don't use this function before calling \ref run(Node) "run()"
379 384
    ///or \ref init(), an instance will be allocated automatically.
380 385
    ///The destructor deallocates this automatically allocated map,
381 386
    ///of course.
382 387
    ///\return <tt> (*this) </tt>
383 388
    Bfs &processedMap(ProcessedMap &m)
384 389
    {
385 390
      if(local_processed) {
386 391
        delete _processed;
387 392
        local_processed=false;
388 393
      }
389 394
      _processed = &m;
390 395
      return *this;
391 396
    }
392 397

	
393 398
    ///Sets the map that stores the distances of the nodes.
394 399

	
395 400
    ///Sets the map that stores the distances of the nodes calculated by
396 401
    ///the algorithm.
397 402
    ///If you don't use this function before calling \ref run(Node) "run()"
398 403
    ///or \ref init(), an instance will be allocated automatically.
399 404
    ///The destructor deallocates this automatically allocated map,
400 405
    ///of course.
401 406
    ///\return <tt> (*this) </tt>
402 407
    Bfs &distMap(DistMap &m)
403 408
    {
404 409
      if(local_dist) {
405 410
        delete _dist;
406 411
        local_dist=false;
407 412
      }
408 413
      _dist = &m;
409 414
      return *this;
410 415
    }
411 416

	
412 417
  public:
413 418

	
414 419
    ///\name Execution Control
415 420
    ///The simplest way to execute the BFS algorithm is to use one of the
416 421
    ///member functions called \ref run(Node) "run()".\n
417 422
    ///If you need better control on the execution, you have to call
418 423
    ///\ref init() first, then you can add several source nodes with
419 424
    ///\ref addSource(). Finally the actual path computation can be
420 425
    ///performed with one of the \ref start() functions.
421 426

	
422 427
    ///@{
423 428

	
424 429
    ///\brief Initializes the internal data structures.
425 430
    ///
426 431
    ///Initializes the internal data structures.
427 432
    void init()
428 433
    {
429 434
      create_maps();
430 435
      _queue.resize(countNodes(*G));
431 436
      _queue_head=_queue_tail=0;
432 437
      _curr_dist=1;
433 438
      for ( NodeIt u(*G) ; u!=INVALID ; ++u ) {
434 439
        _pred->set(u,INVALID);
435 440
        _reached->set(u,false);
436 441
        _processed->set(u,false);
437 442
      }
438 443
    }
439 444

	
440 445
    ///Adds a new source node.
441 446

	
442 447
    ///Adds a new source node to the set of nodes to be processed.
443 448
    ///
444 449
    void addSource(Node s)
445 450
    {
446 451
      if(!(*_reached)[s])
447 452
        {
448 453
          _reached->set(s,true);
449 454
          _pred->set(s,INVALID);
450 455
          _dist->set(s,0);
451 456
          _queue[_queue_head++]=s;
452 457
          _queue_next_dist=_queue_head;
453 458
        }
454 459
    }
455 460

	
456 461
    ///Processes the next node.
457 462

	
458 463
    ///Processes the next node.
459 464
    ///
460 465
    ///\return The processed node.
461 466
    ///
462 467
    ///\pre The queue must not be empty.
463 468
    Node processNextNode()
464 469
    {
465 470
      if(_queue_tail==_queue_next_dist) {
466 471
        _curr_dist++;
467 472
        _queue_next_dist=_queue_head;
468 473
      }
469 474
      Node n=_queue[_queue_tail++];
470 475
      _processed->set(n,true);
471 476
      Node m;
472 477
      for(OutArcIt e(*G,n);e!=INVALID;++e)
473 478
        if(!(*_reached)[m=G->target(e)]) {
474 479
          _queue[_queue_head++]=m;
475 480
          _reached->set(m,true);
476 481
          _pred->set(m,e);
477 482
          _dist->set(m,_curr_dist);
478 483
        }
479 484
      return n;
480 485
    }
481 486

	
482 487
    ///Processes the next node.
483 488

	
484 489
    ///Processes the next node and checks if the given target node
485 490
    ///is reached. If the target node is reachable from the processed
486 491
    ///node, then the \c reach parameter will be set to \c true.
487 492
    ///
488 493
    ///\param target The target node.
489 494
    ///\retval reach Indicates if the target node is reached.
490 495
    ///It should be initially \c false.
491 496
    ///
492 497
    ///\return The processed node.
493 498
    ///
494 499
    ///\pre The queue must not be empty.
495 500
    Node processNextNode(Node target, bool& reach)
496 501
    {
497 502
      if(_queue_tail==_queue_next_dist) {
498 503
        _curr_dist++;
499 504
        _queue_next_dist=_queue_head;
500 505
      }
501 506
      Node n=_queue[_queue_tail++];
502 507
      _processed->set(n,true);
503 508
      Node m;
504 509
      for(OutArcIt e(*G,n);e!=INVALID;++e)
505 510
        if(!(*_reached)[m=G->target(e)]) {
506 511
          _queue[_queue_head++]=m;
507 512
          _reached->set(m,true);
508 513
          _pred->set(m,e);
509 514
          _dist->set(m,_curr_dist);
510 515
          reach = reach || (target == m);
511 516
        }
512 517
      return n;
513 518
    }
514 519

	
515 520
    ///Processes the next node.
516 521

	
517 522
    ///Processes the next node and checks if at least one of reached
518 523
    ///nodes has \c true value in the \c nm node map. If one node
519 524
    ///with \c true value is reachable from the processed node, then the
520 525
    ///\c rnode parameter will be set to the first of such nodes.
521 526
    ///
522 527
    ///\param nm A \c bool (or convertible) node map that indicates the
523 528
    ///possible targets.
524 529
    ///\retval rnode The reached target node.
525 530
    ///It should be initially \c INVALID.
526 531
    ///
527 532
    ///\return The processed node.
528 533
    ///
529 534
    ///\pre The queue must not be empty.
530 535
    template<class NM>
531 536
    Node processNextNode(const NM& nm, Node& rnode)
532 537
    {
533 538
      if(_queue_tail==_queue_next_dist) {
534 539
        _curr_dist++;
535 540
        _queue_next_dist=_queue_head;
536 541
      }
537 542
      Node n=_queue[_queue_tail++];
538 543
      _processed->set(n,true);
539 544
      Node m;
540 545
      for(OutArcIt e(*G,n);e!=INVALID;++e)
541 546
        if(!(*_reached)[m=G->target(e)]) {
542 547
          _queue[_queue_head++]=m;
543 548
          _reached->set(m,true);
544 549
          _pred->set(m,e);
545 550
          _dist->set(m,_curr_dist);
546 551
          if (nm[m] && rnode == INVALID) rnode = m;
547 552
        }
548 553
      return n;
549 554
    }
550 555

	
551 556
    ///The next node to be processed.
552 557

	
553 558
    ///Returns the next node to be processed or \c INVALID if the queue
554 559
    ///is empty.
555 560
    Node nextNode() const
556 561
    {
557 562
      return _queue_tail<_queue_head?_queue[_queue_tail]:INVALID;
558 563
    }
559 564

	
560 565
    ///Returns \c false if there are nodes to be processed.
561 566

	
562 567
    ///Returns \c false if there are nodes to be processed
563 568
    ///in the queue.
564 569
    bool emptyQueue() const { return _queue_tail==_queue_head; }
565 570

	
566 571
    ///Returns the number of the nodes to be processed.
567 572

	
568 573
    ///Returns the number of the nodes to be processed
569 574
    ///in the queue.
570 575
    int queueSize() const { return _queue_head-_queue_tail; }
571 576

	
572 577
    ///Executes the algorithm.
573 578

	
574 579
    ///Executes the algorithm.
575 580
    ///
576 581
    ///This method runs the %BFS algorithm from the root node(s)
577 582
    ///in order to compute the shortest path to each node.
578 583
    ///
579 584
    ///The algorithm computes
580 585
    ///- the shortest path tree (forest),
581 586
    ///- the distance of each node from the root(s).
582 587
    ///
583 588
    ///\pre init() must be called and at least one root node should be
584 589
    ///added with addSource() before using this function.
585 590
    ///
586 591
    ///\note <tt>b.start()</tt> is just a shortcut of the following code.
587 592
    ///\code
588 593
    ///  while ( !b.emptyQueue() ) {
589 594
    ///    b.processNextNode();
590 595
    ///  }
591 596
    ///\endcode
592 597
    void start()
593 598
    {
594 599
      while ( !emptyQueue() ) processNextNode();
595 600
    }
596 601

	
597 602
    ///Executes the algorithm until the given target node is reached.
598 603

	
599 604
    ///Executes the algorithm until the given target node is reached.
600 605
    ///
601 606
    ///This method runs the %BFS algorithm from the root node(s)
602 607
    ///in order to compute the shortest path to \c t.
603 608
    ///
604 609
    ///The algorithm computes
605 610
    ///- the shortest path to \c t,
606 611
    ///- the distance of \c t from the root(s).
607 612
    ///
608 613
    ///\pre init() must be called and at least one root node should be
609 614
    ///added with addSource() before using this function.
610 615
    ///
611 616
    ///\note <tt>b.start(t)</tt> is just a shortcut of the following code.
612 617
    ///\code
613 618
    ///  bool reach = false;
614 619
    ///  while ( !b.emptyQueue() && !reach ) {
615 620
    ///    b.processNextNode(t, reach);
616 621
    ///  }
617 622
    ///\endcode
618 623
    void start(Node t)
619 624
    {
620 625
      bool reach = false;
621 626
      while ( !emptyQueue() && !reach ) processNextNode(t, reach);
622 627
    }
623 628

	
624 629
    ///Executes the algorithm until a condition is met.
625 630

	
626 631
    ///Executes the algorithm until a condition is met.
627 632
    ///
628 633
    ///This method runs the %BFS algorithm from the root node(s) in
629 634
    ///order to compute the shortest path to a node \c v with
630 635
    /// <tt>nm[v]</tt> true, if such a node can be found.
631 636
    ///
632 637
    ///\param nm A \c bool (or convertible) node map. The algorithm
633 638
    ///will stop when it reaches a node \c v with <tt>nm[v]</tt> true.
634 639
    ///
635 640
    ///\return The reached node \c v with <tt>nm[v]</tt> true or
636 641
    ///\c INVALID if no such node was found.
637 642
    ///
638 643
    ///\pre init() must be called and at least one root node should be
639 644
    ///added with addSource() before using this function.
640 645
    ///
641 646
    ///\note <tt>b.start(nm)</tt> is just a shortcut of the following code.
642 647
    ///\code
643 648
    ///  Node rnode = INVALID;
644 649
    ///  while ( !b.emptyQueue() && rnode == INVALID ) {
645 650
    ///    b.processNextNode(nm, rnode);
646 651
    ///  }
647 652
    ///  return rnode;
648 653
    ///\endcode
649 654
    template<class NodeBoolMap>
650 655
    Node start(const NodeBoolMap &nm)
651 656
    {
652 657
      Node rnode = INVALID;
653 658
      while ( !emptyQueue() && rnode == INVALID ) {
654 659
        processNextNode(nm, rnode);
655 660
      }
656 661
      return rnode;
657 662
    }
658 663

	
659 664
    ///Runs the algorithm from the given source node.
660 665

	
661 666
    ///This method runs the %BFS algorithm from node \c s
662 667
    ///in order to compute the shortest path to each node.
663 668
    ///
664 669
    ///The algorithm computes
665 670
    ///- the shortest path tree,
666 671
    ///- the distance of each node from the root.
667 672
    ///
668 673
    ///\note <tt>b.run(s)</tt> is just a shortcut of the following code.
669 674
    ///\code
670 675
    ///  b.init();
671 676
    ///  b.addSource(s);
672 677
    ///  b.start();
673 678
    ///\endcode
674 679
    void run(Node s) {
675 680
      init();
676 681
      addSource(s);
677 682
      start();
678 683
    }
679 684

	
680 685
    ///Finds the shortest path between \c s and \c t.
681 686

	
682 687
    ///This method runs the %BFS algorithm from node \c s
683 688
    ///in order to compute the shortest path to node \c t
684 689
    ///(it stops searching when \c t is processed).
685 690
    ///
686 691
    ///\return \c true if \c t is reachable form \c s.
687 692
    ///
688 693
    ///\note Apart from the return value, <tt>b.run(s,t)</tt> is just a
689 694
    ///shortcut of the following code.
690 695
    ///\code
691 696
    ///  b.init();
692 697
    ///  b.addSource(s);
693 698
    ///  b.start(t);
694 699
    ///\endcode
695 700
    bool run(Node s,Node t) {
696 701
      init();
697 702
      addSource(s);
698 703
      start(t);
699 704
      return reached(t);
700 705
    }
701 706

	
702 707
    ///Runs the algorithm to visit all nodes in the digraph.
703 708

	
704 709
    ///This method runs the %BFS algorithm in order to visit all nodes
705 710
    ///in the digraph.
706 711
    ///
707 712
    ///\note <tt>b.run(s)</tt> is just a shortcut of the following code.
708 713
    ///\code
709 714
    ///  b.init();
710 715
    ///  for (NodeIt n(gr); n != INVALID; ++n) {
711 716
    ///    if (!b.reached(n)) {
712 717
    ///      b.addSource(n);
713 718
    ///      b.start();
714 719
    ///    }
715 720
    ///  }
716 721
    ///\endcode
717 722
    void run() {
718 723
      init();
719 724
      for (NodeIt n(*G); n != INVALID; ++n) {
720 725
        if (!reached(n)) {
721 726
          addSource(n);
722 727
          start();
723 728
        }
724 729
      }
725 730
    }
726 731

	
727 732
    ///@}
728 733

	
729 734
    ///\name Query Functions
730 735
    ///The results of the BFS algorithm can be obtained using these
731 736
    ///functions.\n
732 737
    ///Either \ref run(Node) "run()" or \ref start() should be called
733 738
    ///before using them.
734 739

	
735 740
    ///@{
736 741

	
737 742
    ///The shortest path to the given node.
738 743

	
739 744
    ///Returns the shortest path to the given node from the root(s).
740 745
    ///
741 746
    ///\warning \c t should be reached from the root(s).
742 747
    ///
743 748
    ///\pre Either \ref run(Node) "run()" or \ref init()
744 749
    ///must be called before using this function.
745 750
    Path path(Node t) const { return Path(*G, *_pred, t); }
746 751

	
747 752
    ///The distance of the given node from the root(s).
748 753

	
749 754
    ///Returns the distance of the given node from the root(s).
750 755
    ///
751 756
    ///\warning If node \c v is not reached from the root(s), then
752 757
    ///the return value of this function is undefined.
753 758
    ///
754 759
    ///\pre Either \ref run(Node) "run()" or \ref init()
755 760
    ///must be called before using this function.
756 761
    int dist(Node v) const { return (*_dist)[v]; }
757 762

	
758 763
    ///\brief Returns the 'previous arc' of the shortest path tree for
759 764
    ///the given node.
760 765
    ///
761 766
    ///This function returns the 'previous arc' of the shortest path
762 767
    ///tree for the node \c v, i.e. it returns the last arc of a
763 768
    ///shortest path from a root to \c v. It is \c INVALID if \c v
764 769
    ///is not reached from the root(s) or if \c v is a root.
765 770
    ///
766 771
    ///The shortest path tree used here is equal to the shortest path
767 772
    ///tree used in \ref predNode() and \ref predMap().
768 773
    ///
769 774
    ///\pre Either \ref run(Node) "run()" or \ref init()
770 775
    ///must be called before using this function.
771 776
    Arc predArc(Node v) const { return (*_pred)[v];}
772 777

	
773 778
    ///\brief Returns the 'previous node' of the shortest path tree for
774 779
    ///the given node.
775 780
    ///
776 781
    ///This function returns the 'previous node' of the shortest path
777 782
    ///tree for the node \c v, i.e. it returns the last but one node
778 783
    ///of a shortest path from a root to \c v. It is \c INVALID
779 784
    ///if \c v is not reached from the root(s) or if \c v is a root.
780 785
    ///
781 786
    ///The shortest path tree used here is equal to the shortest path
782 787
    ///tree used in \ref predArc() and \ref predMap().
783 788
    ///
784 789
    ///\pre Either \ref run(Node) "run()" or \ref init()
785 790
    ///must be called before using this function.
786 791
    Node predNode(Node v) const { return (*_pred)[v]==INVALID ? INVALID:
787 792
                                  G->source((*_pred)[v]); }
788 793

	
789 794
    ///\brief Returns a const reference to the node map that stores the
790 795
    /// distances of the nodes.
791 796
    ///
792 797
    ///Returns a const reference to the node map that stores the distances
793 798
    ///of the nodes calculated by the algorithm.
794 799
    ///
795 800
    ///\pre Either \ref run(Node) "run()" or \ref init()
796 801
    ///must be called before using this function.
797 802
    const DistMap &distMap() const { return *_dist;}
798 803

	
799 804
    ///\brief Returns a const reference to the node map that stores the
800 805
    ///predecessor arcs.
801 806
    ///
802 807
    ///Returns a const reference to the node map that stores the predecessor
803 808
    ///arcs, which form the shortest path tree (forest).
804 809
    ///
805 810
    ///\pre Either \ref run(Node) "run()" or \ref init()
806 811
    ///must be called before using this function.
807 812
    const PredMap &predMap() const { return *_pred;}
808 813

	
809 814
    ///Checks if the given node is reached from the root(s).
810 815

	
811 816
    ///Returns \c true if \c v is reached from the root(s).
812 817
    ///
813 818
    ///\pre Either \ref run(Node) "run()" or \ref init()
814 819
    ///must be called before using this function.
815 820
    bool reached(Node v) const { return (*_reached)[v]; }
816 821

	
817 822
    ///@}
818 823
  };
819 824

	
820 825
  ///Default traits class of bfs() function.
821 826

	
822 827
  ///Default traits class of bfs() function.
823 828
  ///\tparam GR Digraph type.
824 829
  template<class GR>
825 830
  struct BfsWizardDefaultTraits
826 831
  {
827 832
    ///The type of the digraph the algorithm runs on.
828 833
    typedef GR Digraph;
829 834

	
830 835
    ///\brief The type of the map that stores the predecessor
831 836
    ///arcs of the shortest paths.
832 837
    ///
833 838
    ///The type of the map that stores the predecessor
834 839
    ///arcs of the shortest paths.
835 840
    ///It must conform to the \ref concepts::WriteMap "WriteMap" concept.
836 841
    typedef typename Digraph::template NodeMap<typename Digraph::Arc> PredMap;
837 842
    ///Instantiates a PredMap.
838 843

	
839 844
    ///This function instantiates a PredMap.
840 845
    ///\param g is the digraph, to which we would like to define the
841 846
    ///PredMap.
842 847
    static PredMap *createPredMap(const Digraph &g)
843 848
    {
844 849
      return new PredMap(g);
845 850
    }
846 851

	
847 852
    ///The type of the map that indicates which nodes are processed.
848 853

	
849 854
    ///The type of the map that indicates which nodes are processed.
850 855
    ///It must conform to the \ref concepts::WriteMap "WriteMap" concept.
851 856
    ///By default, it is a NullMap.
852 857
    typedef NullMap<typename Digraph::Node,bool> ProcessedMap;
853 858
    ///Instantiates a ProcessedMap.
854 859

	
855 860
    ///This function instantiates a ProcessedMap.
856 861
    ///\param g is the digraph, to which
857 862
    ///we would like to define the ProcessedMap.
858 863
#ifdef DOXYGEN
859 864
    static ProcessedMap *createProcessedMap(const Digraph &g)
860 865
#else
861 866
    static ProcessedMap *createProcessedMap(const Digraph &)
862 867
#endif
863 868
    {
864 869
      return new ProcessedMap();
865 870
    }
866 871

	
867 872
    ///The type of the map that indicates which nodes are reached.
868 873

	
869 874
    ///The type of the map that indicates which nodes are reached.
870 875
    ///It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
871 876
    typedef typename Digraph::template NodeMap<bool> ReachedMap;
872 877
    ///Instantiates a ReachedMap.
873 878

	
874 879
    ///This function instantiates a ReachedMap.
875 880
    ///\param g is the digraph, to which
876 881
    ///we would like to define the ReachedMap.
877 882
    static ReachedMap *createReachedMap(const Digraph &g)
878 883
    {
879 884
      return new ReachedMap(g);
880 885
    }
881 886

	
882 887
    ///The type of the map that stores the distances of the nodes.
883 888

	
884 889
    ///The type of the map that stores the distances of the nodes.
885 890
    ///It must conform to the \ref concepts::WriteMap "WriteMap" concept.
886 891
    typedef typename Digraph::template NodeMap<int> DistMap;
887 892
    ///Instantiates a DistMap.
888 893

	
889 894
    ///This function instantiates a DistMap.
890 895
    ///\param g is the digraph, to which we would like to define
891 896
    ///the DistMap
892 897
    static DistMap *createDistMap(const Digraph &g)
893 898
    {
894 899
      return new DistMap(g);
895 900
    }
896 901

	
897 902
    ///The type of the shortest paths.
898 903

	
899 904
    ///The type of the shortest paths.
900 905
    ///It must conform to the \ref concepts::Path "Path" concept.
901 906
    typedef lemon::Path<Digraph> Path;
902 907
  };
903 908

	
904 909
  /// Default traits class used by BfsWizard
905 910

	
906 911
  /// Default traits class used by BfsWizard.
907 912
  /// \tparam GR The type of the digraph.
908 913
  template<class GR>
909 914
  class BfsWizardBase : public BfsWizardDefaultTraits<GR>
910 915
  {
911 916

	
912 917
    typedef BfsWizardDefaultTraits<GR> Base;
913 918
  protected:
914 919
    //The type of the nodes in the digraph.
915 920
    typedef typename Base::Digraph::Node Node;
916 921

	
917 922
    //Pointer to the digraph the algorithm runs on.
918 923
    void *_g;
919 924
    //Pointer to the map of reached nodes.
920 925
    void *_reached;
921 926
    //Pointer to the map of processed nodes.
922 927
    void *_processed;
923 928
    //Pointer to the map of predecessors arcs.
924 929
    void *_pred;
925 930
    //Pointer to the map of distances.
926 931
    void *_dist;
927 932
    //Pointer to the shortest path to the target node.
928 933
    void *_path;
929 934
    //Pointer to the distance of the target node.
930 935
    int *_di;
931 936

	
932 937
    public:
933 938
    /// Constructor.
934 939

	
935 940
    /// This constructor does not require parameters, it initiates
936 941
    /// all of the attributes to \c 0.
937 942
    BfsWizardBase() : _g(0), _reached(0), _processed(0), _pred(0),
938 943
                      _dist(0), _path(0), _di(0) {}
939 944

	
940 945
    /// Constructor.
941 946

	
942 947
    /// This constructor requires one parameter,
943 948
    /// others are initiated to \c 0.
944 949
    /// \param g The digraph the algorithm runs on.
945 950
    BfsWizardBase(const GR &g) :
946 951
      _g(reinterpret_cast<void*>(const_cast<GR*>(&g))),
947 952
      _reached(0), _processed(0), _pred(0), _dist(0),  _path(0), _di(0) {}
948 953

	
949 954
  };
950 955

	
951 956
  /// Auxiliary class for the function-type interface of BFS algorithm.
952 957

	
953 958
  /// This auxiliary class is created to implement the
954 959
  /// \ref bfs() "function-type interface" of \ref Bfs algorithm.
955 960
  /// It does not have own \ref run(Node) "run()" method, it uses the
956 961
  /// functions and features of the plain \ref Bfs.
957 962
  ///
958 963
  /// This class should only be used through the \ref bfs() function,
959 964
  /// which makes it easier to use the algorithm.
965
  ///
966
  /// \tparam TR The traits class that defines various types used by the
967
  /// algorithm.
960 968
  template<class TR>
961 969
  class BfsWizard : public TR
962 970
  {
963 971
    typedef TR Base;
964 972

	
965 973
    typedef typename TR::Digraph Digraph;
966 974

	
967 975
    typedef typename Digraph::Node Node;
968 976
    typedef typename Digraph::NodeIt NodeIt;
969 977
    typedef typename Digraph::Arc Arc;
970 978
    typedef typename Digraph::OutArcIt OutArcIt;
971 979

	
972 980
    typedef typename TR::PredMap PredMap;
973 981
    typedef typename TR::DistMap DistMap;
974 982
    typedef typename TR::ReachedMap ReachedMap;
975 983
    typedef typename TR::ProcessedMap ProcessedMap;
976 984
    typedef typename TR::Path Path;
977 985

	
978 986
  public:
979 987

	
980 988
    /// Constructor.
981 989
    BfsWizard() : TR() {}
982 990

	
983 991
    /// Constructor that requires parameters.
984 992

	
985 993
    /// Constructor that requires parameters.
986 994
    /// These parameters will be the default values for the traits class.
987 995
    /// \param g The digraph the algorithm runs on.
988 996
    BfsWizard(const Digraph &g) :
989 997
      TR(g) {}
990 998

	
991 999
    ///Copy constructor
992 1000
    BfsWizard(const TR &b) : TR(b) {}
993 1001

	
994 1002
    ~BfsWizard() {}
995 1003

	
996 1004
    ///Runs BFS algorithm from the given source node.
997 1005

	
998 1006
    ///This method runs BFS algorithm from node \c s
999 1007
    ///in order to compute the shortest path to each node.
1000 1008
    void run(Node s)
1001 1009
    {
1002 1010
      Bfs<Digraph,TR> alg(*reinterpret_cast<const Digraph*>(Base::_g));
1003 1011
      if (Base::_pred)
1004 1012
        alg.predMap(*reinterpret_cast<PredMap*>(Base::_pred));
1005 1013
      if (Base::_dist)
1006 1014
        alg.distMap(*reinterpret_cast<DistMap*>(Base::_dist));
1007 1015
      if (Base::_reached)
1008 1016
        alg.reachedMap(*reinterpret_cast<ReachedMap*>(Base::_reached));
1009 1017
      if (Base::_processed)
1010 1018
        alg.processedMap(*reinterpret_cast<ProcessedMap*>(Base::_processed));
1011 1019
      if (s!=INVALID)
1012 1020
        alg.run(s);
1013 1021
      else
1014 1022
        alg.run();
1015 1023
    }
1016 1024

	
1017 1025
    ///Finds the shortest path between \c s and \c t.
1018 1026

	
1019 1027
    ///This method runs BFS algorithm from node \c s
1020 1028
    ///in order to compute the shortest path to node \c t
1021 1029
    ///(it stops searching when \c t is processed).
1022 1030
    ///
1023 1031
    ///\return \c true if \c t is reachable form \c s.
1024 1032
    bool run(Node s, Node t)
1025 1033
    {
1026 1034
      Bfs<Digraph,TR> alg(*reinterpret_cast<const Digraph*>(Base::_g));
1027 1035
      if (Base::_pred)
1028 1036
        alg.predMap(*reinterpret_cast<PredMap*>(Base::_pred));
1029 1037
      if (Base::_dist)
1030 1038
        alg.distMap(*reinterpret_cast<DistMap*>(Base::_dist));
1031 1039
      if (Base::_reached)
1032 1040
        alg.reachedMap(*reinterpret_cast<ReachedMap*>(Base::_reached));
1033 1041
      if (Base::_processed)
1034 1042
        alg.processedMap(*reinterpret_cast<ProcessedMap*>(Base::_processed));
1035 1043
      alg.run(s,t);
1036 1044
      if (Base::_path)
1037 1045
        *reinterpret_cast<Path*>(Base::_path) = alg.path(t);
1038 1046
      if (Base::_di)
1039 1047
        *Base::_di = alg.dist(t);
1040 1048
      return alg.reached(t);
1041 1049
    }
1042 1050

	
1043 1051
    ///Runs BFS algorithm to visit all nodes in the digraph.
1044 1052

	
1045 1053
    ///This method runs BFS algorithm in order to visit all nodes
1046 1054
    ///in the digraph.
1047 1055
    void run()
1048 1056
    {
1049 1057
      run(INVALID);
1050 1058
    }
1051 1059

	
1052 1060
    template<class T>
1053 1061
    struct SetPredMapBase : public Base {
1054 1062
      typedef T PredMap;
1055 1063
      static PredMap *createPredMap(const Digraph &) { return 0; };
1056 1064
      SetPredMapBase(const TR &b) : TR(b) {}
1057 1065
    };
1058 1066

	
1059 1067
    ///\brief \ref named-templ-param "Named parameter" for setting
1060 1068
    ///the predecessor map.
1061 1069
    ///
1062 1070
    ///\ref named-templ-param "Named parameter" function for setting
1063 1071
    ///the map that stores the predecessor arcs of the nodes.
1064 1072
    template<class T>
1065 1073
    BfsWizard<SetPredMapBase<T> > predMap(const T &t)
1066 1074
    {
1067 1075
      Base::_pred=reinterpret_cast<void*>(const_cast<T*>(&t));
1068 1076
      return BfsWizard<SetPredMapBase<T> >(*this);
1069 1077
    }
1070 1078

	
1071 1079
    template<class T>
1072 1080
    struct SetReachedMapBase : public Base {
1073 1081
      typedef T ReachedMap;
1074 1082
      static ReachedMap *createReachedMap(const Digraph &) { return 0; };
1075 1083
      SetReachedMapBase(const TR &b) : TR(b) {}
1076 1084
    };
1077 1085

	
1078 1086
    ///\brief \ref named-templ-param "Named parameter" for setting
1079 1087
    ///the reached map.
1080 1088
    ///
1081 1089
    ///\ref named-templ-param "Named parameter" function for setting
1082 1090
    ///the map that indicates which nodes are reached.
1083 1091
    template<class T>
1084 1092
    BfsWizard<SetReachedMapBase<T> > reachedMap(const T &t)
1085 1093
    {
1086 1094
      Base::_reached=reinterpret_cast<void*>(const_cast<T*>(&t));
1087 1095
      return BfsWizard<SetReachedMapBase<T> >(*this);
1088 1096
    }
1089 1097

	
1090 1098
    template<class T>
1091 1099
    struct SetDistMapBase : public Base {
1092 1100
      typedef T DistMap;
1093 1101
      static DistMap *createDistMap(const Digraph &) { return 0; };
1094 1102
      SetDistMapBase(const TR &b) : TR(b) {}
1095 1103
    };
1096 1104

	
1097 1105
    ///\brief \ref named-templ-param "Named parameter" for setting
1098 1106
    ///the distance map.
1099 1107
    ///
1100 1108
    ///\ref named-templ-param "Named parameter" function for setting
1101 1109
    ///the map that stores the distances of the nodes calculated
1102 1110
    ///by the algorithm.
1103 1111
    template<class T>
1104 1112
    BfsWizard<SetDistMapBase<T> > distMap(const T &t)
1105 1113
    {
1106 1114
      Base::_dist=reinterpret_cast<void*>(const_cast<T*>(&t));
1107 1115
      return BfsWizard<SetDistMapBase<T> >(*this);
1108 1116
    }
1109 1117

	
1110 1118
    template<class T>
1111 1119
    struct SetProcessedMapBase : public Base {
1112 1120
      typedef T ProcessedMap;
1113 1121
      static ProcessedMap *createProcessedMap(const Digraph &) { return 0; };
1114 1122
      SetProcessedMapBase(const TR &b) : TR(b) {}
1115 1123
    };
1116 1124

	
1117 1125
    ///\brief \ref named-func-param "Named parameter" for setting
1118 1126
    ///the processed map.
1119 1127
    ///
1120 1128
    ///\ref named-templ-param "Named parameter" function for setting
1121 1129
    ///the map that indicates which nodes are processed.
1122 1130
    template<class T>
1123 1131
    BfsWizard<SetProcessedMapBase<T> > processedMap(const T &t)
1124 1132
    {
1125 1133
      Base::_processed=reinterpret_cast<void*>(const_cast<T*>(&t));
1126 1134
      return BfsWizard<SetProcessedMapBase<T> >(*this);
1127 1135
    }
1128 1136

	
1129 1137
    template<class T>
1130 1138
    struct SetPathBase : public Base {
1131 1139
      typedef T Path;
1132 1140
      SetPathBase(const TR &b) : TR(b) {}
1133 1141
    };
1134 1142
    ///\brief \ref named-func-param "Named parameter"
1135 1143
    ///for getting the shortest path to the target node.
1136 1144
    ///
1137 1145
    ///\ref named-func-param "Named parameter"
1138 1146
    ///for getting the shortest path to the target node.
1139 1147
    template<class T>
1140 1148
    BfsWizard<SetPathBase<T> > path(const T &t)
1141 1149
    {
1142 1150
      Base::_path=reinterpret_cast<void*>(const_cast<T*>(&t));
1143 1151
      return BfsWizard<SetPathBase<T> >(*this);
1144 1152
    }
1145 1153

	
1146 1154
    ///\brief \ref named-func-param "Named parameter"
1147 1155
    ///for getting the distance of the target node.
1148 1156
    ///
1149 1157
    ///\ref named-func-param "Named parameter"
1150 1158
    ///for getting the distance of the target node.
1151 1159
    BfsWizard dist(const int &d)
1152 1160
    {
1153 1161
      Base::_di=const_cast<int*>(&d);
1154 1162
      return *this;
1155 1163
    }
1156 1164

	
1157 1165
  };
1158 1166

	
1159 1167
  ///Function-type interface for BFS algorithm.
1160 1168

	
1161 1169
  /// \ingroup search
1162 1170
  ///Function-type interface for BFS algorithm.
1163 1171
  ///
1164 1172
  ///This function also has several \ref named-func-param "named parameters",
1165 1173
  ///they are declared as the members of class \ref BfsWizard.
1166 1174
  ///The following examples show how to use these parameters.
1167 1175
  ///\code
1168 1176
  ///  // Compute shortest path from node s to each node
1169 1177
  ///  bfs(g).predMap(preds).distMap(dists).run(s);
1170 1178
  ///
1171 1179
  ///  // Compute shortest path from s to t
1172 1180
  ///  bool reached = bfs(g).path(p).dist(d).run(s,t);
1173 1181
  ///\endcode
1174 1182
  ///\warning Don't forget to put the \ref BfsWizard::run(Node) "run()"
1175 1183
  ///to the end of the parameter list.
1176 1184
  ///\sa BfsWizard
1177 1185
  ///\sa Bfs
1178 1186
  template<class GR>
1179 1187
  BfsWizard<BfsWizardBase<GR> >
1180 1188
  bfs(const GR &digraph)
1181 1189
  {
1182 1190
    return BfsWizard<BfsWizardBase<GR> >(digraph);
1183 1191
  }
1184 1192

	
1185 1193
#ifdef DOXYGEN
1186 1194
  /// \brief Visitor class for BFS.
1187 1195
  ///
1188 1196
  /// This class defines the interface of the BfsVisit events, and
1189 1197
  /// it could be the base of a real visitor class.
1190 1198
  template <typename GR>
1191 1199
  struct BfsVisitor {
1192 1200
    typedef GR Digraph;
1193 1201
    typedef typename Digraph::Arc Arc;
1194 1202
    typedef typename Digraph::Node Node;
1195 1203
    /// \brief Called for the source node(s) of the BFS.
1196 1204
    ///
1197 1205
    /// This function is called for the source node(s) of the BFS.
1198 1206
    void start(const Node& node) {}
1199 1207
    /// \brief Called when a node is reached first time.
1200 1208
    ///
1201 1209
    /// This function is called when a node is reached first time.
1202 1210
    void reach(const Node& node) {}
1203 1211
    /// \brief Called when a node is processed.
1204 1212
    ///
1205 1213
    /// This function is called when a node is processed.
1206 1214
    void process(const Node& node) {}
1207 1215
    /// \brief Called when an arc reaches a new node.
1208 1216
    ///
1209 1217
    /// This function is called when the BFS finds an arc whose target node
1210 1218
    /// is not reached yet.
1211 1219
    void discover(const Arc& arc) {}
1212 1220
    /// \brief Called when an arc is examined but its target node is
1213 1221
    /// already discovered.
1214 1222
    ///
1215 1223
    /// This function is called when an arc is examined but its target node is
1216 1224
    /// already discovered.
1217 1225
    void examine(const Arc& arc) {}
1218 1226
  };
1219 1227
#else
1220 1228
  template <typename GR>
1221 1229
  struct BfsVisitor {
1222 1230
    typedef GR Digraph;
1223 1231
    typedef typename Digraph::Arc Arc;
1224 1232
    typedef typename Digraph::Node Node;
1225 1233
    void start(const Node&) {}
1226 1234
    void reach(const Node&) {}
1227 1235
    void process(const Node&) {}
1228 1236
    void discover(const Arc&) {}
1229 1237
    void examine(const Arc&) {}
1230 1238

	
1231 1239
    template <typename _Visitor>
1232 1240
    struct Constraints {
1233 1241
      void constraints() {
1234 1242
        Arc arc;
1235 1243
        Node node;
1236 1244
        visitor.start(node);
1237 1245
        visitor.reach(node);
1238 1246
        visitor.process(node);
1239 1247
        visitor.discover(arc);
1240 1248
        visitor.examine(arc);
1241 1249
      }
1242 1250
      _Visitor& visitor;
1243 1251
    };
1244 1252
  };
1245 1253
#endif
1246 1254

	
1247 1255
  /// \brief Default traits class of BfsVisit class.
1248 1256
  ///
1249 1257
  /// Default traits class of BfsVisit class.
1250 1258
  /// \tparam GR The type of the digraph the algorithm runs on.
1251 1259
  template<class GR>
1252 1260
  struct BfsVisitDefaultTraits {
1253 1261

	
1254 1262
    /// \brief The type of the digraph the algorithm runs on.
1255 1263
    typedef GR Digraph;
1256 1264

	
1257 1265
    /// \brief The type of the map that indicates which nodes are reached.
1258 1266
    ///
1259 1267
    /// The type of the map that indicates which nodes are reached.
1260 1268
    /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
1261 1269
    typedef typename Digraph::template NodeMap<bool> ReachedMap;
1262 1270

	
1263 1271
    /// \brief Instantiates a ReachedMap.
1264 1272
    ///
1265 1273
    /// This function instantiates a ReachedMap.
1266 1274
    /// \param digraph is the digraph, to which
1267 1275
    /// we would like to define the ReachedMap.
1268 1276
    static ReachedMap *createReachedMap(const Digraph &digraph) {
1269 1277
      return new ReachedMap(digraph);
1270 1278
    }
1271 1279

	
1272 1280
  };
1273 1281

	
1274 1282
  /// \ingroup search
1275 1283
  ///
1276 1284
  /// \brief BFS algorithm class with visitor interface.
1277 1285
  ///
1278 1286
  /// This class provides an efficient implementation of the BFS algorithm
1279 1287
  /// with visitor interface.
1280 1288
  ///
1281 1289
  /// The BfsVisit class provides an alternative interface to the Bfs
1282 1290
  /// class. It works with callback mechanism, the BfsVisit object calls
1283 1291
  /// the member functions of the \c Visitor class on every BFS event.
1284 1292
  ///
1285 1293
  /// This interface of the BFS algorithm should be used in special cases
1286 1294
  /// when extra actions have to be performed in connection with certain
1287 1295
  /// events of the BFS algorithm. Otherwise consider to use Bfs or bfs()
1288 1296
  /// instead.
1289 1297
  ///
1290 1298
  /// \tparam GR The type of the digraph the algorithm runs on.
1291 1299
  /// The default type is \ref ListDigraph.
1292 1300
  /// The value of GR is not used directly by \ref BfsVisit,
1293 1301
  /// it is only passed to \ref BfsVisitDefaultTraits.
1294 1302
  /// \tparam VS The Visitor type that is used by the algorithm.
1295 1303
  /// \ref BfsVisitor "BfsVisitor<GR>" is an empty visitor, which
1296 1304
  /// does not observe the BFS events. If you want to observe the BFS
1297 1305
  /// events, you should implement your own visitor class.
1298
  /// \tparam TR Traits class to set various data types used by the
1299
  /// algorithm. The default traits class is
1300
  /// \ref BfsVisitDefaultTraits "BfsVisitDefaultTraits<GR>".
1301
  /// See \ref BfsVisitDefaultTraits for the documentation of
1302
  /// a BFS visit traits class.
1306
  /// \tparam TR The traits class that defines various types used by the
1307
  /// algorithm. By default, it is \ref BfsVisitDefaultTraits
1308
  /// "BfsVisitDefaultTraits<GR>".
1309
  /// In most cases, this parameter should not be set directly,
1310
  /// consider to use the named template parameters instead.
1303 1311
#ifdef DOXYGEN
1304 1312
  template <typename GR, typename VS, typename TR>
1305 1313
#else
1306 1314
  template <typename GR = ListDigraph,
1307 1315
            typename VS = BfsVisitor<GR>,
1308 1316
            typename TR = BfsVisitDefaultTraits<GR> >
1309 1317
#endif
1310 1318
  class BfsVisit {
1311 1319
  public:
1312 1320

	
1313 1321
    ///The traits class.
1314 1322
    typedef TR Traits;
1315 1323

	
1316 1324
    ///The type of the digraph the algorithm runs on.
1317 1325
    typedef typename Traits::Digraph Digraph;
1318 1326

	
1319 1327
    ///The visitor type used by the algorithm.
1320 1328
    typedef VS Visitor;
1321 1329

	
1322 1330
    ///The type of the map that indicates which nodes are reached.
1323 1331
    typedef typename Traits::ReachedMap ReachedMap;
1324 1332

	
1325 1333
  private:
1326 1334

	
1327 1335
    typedef typename Digraph::Node Node;
1328 1336
    typedef typename Digraph::NodeIt NodeIt;
1329 1337
    typedef typename Digraph::Arc Arc;
1330 1338
    typedef typename Digraph::OutArcIt OutArcIt;
1331 1339

	
1332 1340
    //Pointer to the underlying digraph.
1333 1341
    const Digraph *_digraph;
1334 1342
    //Pointer to the visitor object.
1335 1343
    Visitor *_visitor;
1336 1344
    //Pointer to the map of reached status of the nodes.
1337 1345
    ReachedMap *_reached;
1338 1346
    //Indicates if _reached is locally allocated (true) or not.
1339 1347
    bool local_reached;
1340 1348

	
1341 1349
    std::vector<typename Digraph::Node> _list;
1342 1350
    int _list_front, _list_back;
1343 1351

	
1344 1352
    //Creates the maps if necessary.
1345 1353
    void create_maps() {
1346 1354
      if(!_reached) {
1347 1355
        local_reached = true;
1348 1356
        _reached = Traits::createReachedMap(*_digraph);
1349 1357
      }
1350 1358
    }
1351 1359

	
1352 1360
  protected:
1353 1361

	
1354 1362
    BfsVisit() {}
1355 1363

	
1356 1364
  public:
1357 1365

	
1358 1366
    typedef BfsVisit Create;
1359 1367

	
1360 1368
    /// \name Named Template Parameters
1361 1369

	
1362 1370
    ///@{
1363 1371
    template <class T>
1364 1372
    struct SetReachedMapTraits : public Traits {
1365 1373
      typedef T ReachedMap;
1366 1374
      static ReachedMap *createReachedMap(const Digraph &digraph) {
1367 1375
        LEMON_ASSERT(false, "ReachedMap is not initialized");
1368 1376
        return 0; // ignore warnings
1369 1377
      }
1370 1378
    };
1371 1379
    /// \brief \ref named-templ-param "Named parameter" for setting
1372 1380
    /// ReachedMap type.
1373 1381
    ///
1374 1382
    /// \ref named-templ-param "Named parameter" for setting ReachedMap type.
1375 1383
    template <class T>
1376 1384
    struct SetReachedMap : public BfsVisit< Digraph, Visitor,
1377 1385
                                            SetReachedMapTraits<T> > {
1378 1386
      typedef BfsVisit< Digraph, Visitor, SetReachedMapTraits<T> > Create;
1379 1387
    };
1380 1388
    ///@}
1381 1389

	
1382 1390
  public:
1383 1391

	
1384 1392
    /// \brief Constructor.
1385 1393
    ///
1386 1394
    /// Constructor.
1387 1395
    ///
1388 1396
    /// \param digraph The digraph the algorithm runs on.
1389 1397
    /// \param visitor The visitor object of the algorithm.
1390 1398
    BfsVisit(const Digraph& digraph, Visitor& visitor)
1391 1399
      : _digraph(&digraph), _visitor(&visitor),
1392 1400
        _reached(0), local_reached(false) {}
1393 1401

	
1394 1402
    /// \brief Destructor.
1395 1403
    ~BfsVisit() {
1396 1404
      if(local_reached) delete _reached;
1397 1405
    }
1398 1406

	
1399 1407
    /// \brief Sets the map that indicates which nodes are reached.
1400 1408
    ///
1401 1409
    /// Sets the map that indicates which nodes are reached.
1402 1410
    /// If you don't use this function before calling \ref run(Node) "run()"
1403 1411
    /// or \ref init(), an instance will be allocated automatically.
1404 1412
    /// The destructor deallocates this automatically allocated map,
1405 1413
    /// of course.
1406 1414
    /// \return <tt> (*this) </tt>
1407 1415
    BfsVisit &reachedMap(ReachedMap &m) {
1408 1416
      if(local_reached) {
1409 1417
        delete _reached;
1410 1418
        local_reached = false;
1411 1419
      }
1412 1420
      _reached = &m;
1413 1421
      return *this;
1414 1422
    }
1415 1423

	
1416 1424
  public:
1417 1425

	
1418 1426
    /// \name Execution Control
1419 1427
    /// The simplest way to execute the BFS algorithm is to use one of the
1420 1428
    /// member functions called \ref run(Node) "run()".\n
1421 1429
    /// If you need better control on the execution, you have to call
1422 1430
    /// \ref init() first, then you can add several source nodes with
1423 1431
    /// \ref addSource(). Finally the actual path computation can be
1424 1432
    /// performed with one of the \ref start() functions.
1425 1433

	
1426 1434
    /// @{
1427 1435

	
1428 1436
    /// \brief Initializes the internal data structures.
1429 1437
    ///
1430 1438
    /// Initializes the internal data structures.
1431 1439
    void init() {
1432 1440
      create_maps();
1433 1441
      _list.resize(countNodes(*_digraph));
1434 1442
      _list_front = _list_back = -1;
1435 1443
      for (NodeIt u(*_digraph) ; u != INVALID ; ++u) {
1436 1444
        _reached->set(u, false);
1437 1445
      }
1438 1446
    }
1439 1447

	
1440 1448
    /// \brief Adds a new source node.
1441 1449
    ///
1442 1450
    /// Adds a new source node to the set of nodes to be processed.
1443 1451
    void addSource(Node s) {
1444 1452
      if(!(*_reached)[s]) {
1445 1453
          _reached->set(s,true);
1446 1454
          _visitor->start(s);
1447 1455
          _visitor->reach(s);
1448 1456
          _list[++_list_back] = s;
1449 1457
        }
1450 1458
    }
1451 1459

	
1452 1460
    /// \brief Processes the next node.
1453 1461
    ///
1454 1462
    /// Processes the next node.
1455 1463
    ///
1456 1464
    /// \return The processed node.
1457 1465
    ///
1458 1466
    /// \pre The queue must not be empty.
1459 1467
    Node processNextNode() {
1460 1468
      Node n = _list[++_list_front];
1461 1469
      _visitor->process(n);
1462 1470
      Arc e;
1463 1471
      for (_digraph->firstOut(e, n); e != INVALID; _digraph->nextOut(e)) {
1464 1472
        Node m = _digraph->target(e);
1465 1473
        if (!(*_reached)[m]) {
1466 1474
          _visitor->discover(e);
1467 1475
          _visitor->reach(m);
1468 1476
          _reached->set(m, true);
1469 1477
          _list[++_list_back] = m;
1470 1478
        } else {
1471 1479
          _visitor->examine(e);
1472 1480
        }
1473 1481
      }
1474 1482
      return n;
1475 1483
    }
1476 1484

	
1477 1485
    /// \brief Processes the next node.
1478 1486
    ///
1479 1487
    /// Processes the next node and checks if the given target node
1480 1488
    /// is reached. If the target node is reachable from the processed
1481 1489
    /// node, then the \c reach parameter will be set to \c true.
1482 1490
    ///
1483 1491
    /// \param target The target node.
1484 1492
    /// \retval reach Indicates if the target node is reached.
1485 1493
    /// It should be initially \c false.
1486 1494
    ///
1487 1495
    /// \return The processed node.
1488 1496
    ///
1489 1497
    /// \pre The queue must not be empty.
1490 1498
    Node processNextNode(Node target, bool& reach) {
1491 1499
      Node n = _list[++_list_front];
1492 1500
      _visitor->process(n);
1493 1501
      Arc e;
1494 1502
      for (_digraph->firstOut(e, n); e != INVALID; _digraph->nextOut(e)) {
1495 1503
        Node m = _digraph->target(e);
1496 1504
        if (!(*_reached)[m]) {
1497 1505
          _visitor->discover(e);
1498 1506
          _visitor->reach(m);
1499 1507
          _reached->set(m, true);
1500 1508
          _list[++_list_back] = m;
1501 1509
          reach = reach || (target == m);
1502 1510
        } else {
1503 1511
          _visitor->examine(e);
1504 1512
        }
1505 1513
      }
1506 1514
      return n;
1507 1515
    }
1508 1516

	
1509 1517
    /// \brief Processes the next node.
1510 1518
    ///
1511 1519
    /// Processes the next node and checks if at least one of reached
1512 1520
    /// nodes has \c true value in the \c nm node map. If one node
1513 1521
    /// with \c true value is reachable from the processed node, then the
1514 1522
    /// \c rnode parameter will be set to the first of such nodes.
1515 1523
    ///
1516 1524
    /// \param nm A \c bool (or convertible) node map that indicates the
1517 1525
    /// possible targets.
1518 1526
    /// \retval rnode The reached target node.
1519 1527
    /// It should be initially \c INVALID.
1520 1528
    ///
1521 1529
    /// \return The processed node.
1522 1530
    ///
1523 1531
    /// \pre The queue must not be empty.
1524 1532
    template <typename NM>
1525 1533
    Node processNextNode(const NM& nm, Node& rnode) {
1526 1534
      Node n = _list[++_list_front];
1527 1535
      _visitor->process(n);
1528 1536
      Arc e;
1529 1537
      for (_digraph->firstOut(e, n); e != INVALID; _digraph->nextOut(e)) {
1530 1538
        Node m = _digraph->target(e);
1531 1539
        if (!(*_reached)[m]) {
1532 1540
          _visitor->discover(e);
1533 1541
          _visitor->reach(m);
1534 1542
          _reached->set(m, true);
1535 1543
          _list[++_list_back] = m;
1536 1544
          if (nm[m] && rnode == INVALID) rnode = m;
1537 1545
        } else {
1538 1546
          _visitor->examine(e);
1539 1547
        }
1540 1548
      }
1541 1549
      return n;
1542 1550
    }
1543 1551

	
1544 1552
    /// \brief The next node to be processed.
1545 1553
    ///
1546 1554
    /// Returns the next node to be processed or \c INVALID if the queue
1547 1555
    /// is empty.
1548 1556
    Node nextNode() const {
1549 1557
      return _list_front != _list_back ? _list[_list_front + 1] : INVALID;
1550 1558
    }
1551 1559

	
1552 1560
    /// \brief Returns \c false if there are nodes
1553 1561
    /// to be processed.
1554 1562
    ///
1555 1563
    /// Returns \c false if there are nodes
1556 1564
    /// to be processed in the queue.
1557 1565
    bool emptyQueue() const { return _list_front == _list_back; }
1558 1566

	
1559 1567
    /// \brief Returns the number of the nodes to be processed.
1560 1568
    ///
1561 1569
    /// Returns the number of the nodes to be processed in the queue.
1562 1570
    int queueSize() const { return _list_back - _list_front; }
1563 1571

	
1564 1572
    /// \brief Executes the algorithm.
1565 1573
    ///
1566 1574
    /// Executes the algorithm.
1567 1575
    ///
1568 1576
    /// This method runs the %BFS algorithm from the root node(s)
1569 1577
    /// in order to compute the shortest path to each node.
1570 1578
    ///
1571 1579
    /// The algorithm computes
1572 1580
    /// - the shortest path tree (forest),
1573 1581
    /// - the distance of each node from the root(s).
1574 1582
    ///
1575 1583
    /// \pre init() must be called and at least one root node should be added
1576 1584
    /// with addSource() before using this function.
1577 1585
    ///
1578 1586
    /// \note <tt>b.start()</tt> is just a shortcut of the following code.
1579 1587
    /// \code
1580 1588
    ///   while ( !b.emptyQueue() ) {
1581 1589
    ///     b.processNextNode();
1582 1590
    ///   }
1583 1591
    /// \endcode
1584 1592
    void start() {
1585 1593
      while ( !emptyQueue() ) processNextNode();
1586 1594
    }
1587 1595

	
1588 1596
    /// \brief Executes the algorithm until the given target node is reached.
1589 1597
    ///
1590 1598
    /// Executes the algorithm until the given target node is reached.
1591 1599
    ///
1592 1600
    /// This method runs the %BFS algorithm from the root node(s)
1593 1601
    /// in order to compute the shortest path to \c t.
1594 1602
    ///
1595 1603
    /// The algorithm computes
1596 1604
    /// - the shortest path to \c t,
1597 1605
    /// - the distance of \c t from the root(s).
1598 1606
    ///
1599 1607
    /// \pre init() must be called and at least one root node should be
1600 1608
    /// added with addSource() before using this function.
1601 1609
    ///
1602 1610
    /// \note <tt>b.start(t)</tt> is just a shortcut of the following code.
1603 1611
    /// \code
1604 1612
    ///   bool reach = false;
1605 1613
    ///   while ( !b.emptyQueue() && !reach ) {
1606 1614
    ///     b.processNextNode(t, reach);
1607 1615
    ///   }
1608 1616
    /// \endcode
1609 1617
    void start(Node t) {
1610 1618
      bool reach = false;
1611 1619
      while ( !emptyQueue() && !reach ) processNextNode(t, reach);
1612 1620
    }
1613 1621

	
1614 1622
    /// \brief Executes the algorithm until a condition is met.
1615 1623
    ///
1616 1624
    /// Executes the algorithm until a condition is met.
1617 1625
    ///
1618 1626
    /// This method runs the %BFS algorithm from the root node(s) in
1619 1627
    /// order to compute the shortest path to a node \c v with
1620 1628
    /// <tt>nm[v]</tt> true, if such a node can be found.
1621 1629
    ///
1622 1630
    /// \param nm must be a bool (or convertible) node map. The
1623 1631
    /// algorithm will stop when it reaches a node \c v with
1624 1632
    /// <tt>nm[v]</tt> true.
1625 1633
    ///
1626 1634
    /// \return The reached node \c v with <tt>nm[v]</tt> true or
1627 1635
    /// \c INVALID if no such node was found.
1628 1636
    ///
1629 1637
    /// \pre init() must be called and at least one root node should be
1630 1638
    /// added with addSource() before using this function.
1631 1639
    ///
1632 1640
    /// \note <tt>b.start(nm)</tt> is just a shortcut of the following code.
1633 1641
    /// \code
1634 1642
    ///   Node rnode = INVALID;
1635 1643
    ///   while ( !b.emptyQueue() && rnode == INVALID ) {
1636 1644
    ///     b.processNextNode(nm, rnode);
1637 1645
    ///   }
1638 1646
    ///   return rnode;
1639 1647
    /// \endcode
1640 1648
    template <typename NM>
1641 1649
    Node start(const NM &nm) {
1642 1650
      Node rnode = INVALID;
1643 1651
      while ( !emptyQueue() && rnode == INVALID ) {
1644 1652
        processNextNode(nm, rnode);
1645 1653
      }
1646 1654
      return rnode;
1647 1655
    }
1648 1656

	
1649 1657
    /// \brief Runs the algorithm from the given source node.
1650 1658
    ///
1651 1659
    /// This method runs the %BFS algorithm from node \c s
1652 1660
    /// in order to compute the shortest path to each node.
1653 1661
    ///
1654 1662
    /// The algorithm computes
1655 1663
    /// - the shortest path tree,
1656 1664
    /// - the distance of each node from the root.
1657 1665
    ///
1658 1666
    /// \note <tt>b.run(s)</tt> is just a shortcut of the following code.
1659 1667
    ///\code
1660 1668
    ///   b.init();
1661 1669
    ///   b.addSource(s);
1662 1670
    ///   b.start();
1663 1671
    ///\endcode
1664 1672
    void run(Node s) {
1665 1673
      init();
1666 1674
      addSource(s);
1667 1675
      start();
1668 1676
    }
1669 1677

	
1670 1678
    /// \brief Finds the shortest path between \c s and \c t.
1671 1679
    ///
1672 1680
    /// This method runs the %BFS algorithm from node \c s
1673 1681
    /// in order to compute the shortest path to node \c t
1674 1682
    /// (it stops searching when \c t is processed).
1675 1683
    ///
1676 1684
    /// \return \c true if \c t is reachable form \c s.
1677 1685
    ///
1678 1686
    /// \note Apart from the return value, <tt>b.run(s,t)</tt> is just a
1679 1687
    /// shortcut of the following code.
1680 1688
    ///\code
1681 1689
    ///   b.init();
1682 1690
    ///   b.addSource(s);
1683 1691
    ///   b.start(t);
1684 1692
    ///\endcode
1685 1693
    bool run(Node s,Node t) {
1686 1694
      init();
1687 1695
      addSource(s);
1688 1696
      start(t);
1689 1697
      return reached(t);
1690 1698
    }
1691 1699

	
1692 1700
    /// \brief Runs the algorithm to visit all nodes in the digraph.
1693 1701
    ///
1694 1702
    /// This method runs the %BFS algorithm in order to visit all nodes
1695 1703
    /// in the digraph.
1696 1704
    ///
1697 1705
    /// \note <tt>b.run(s)</tt> is just a shortcut of the following code.
1698 1706
    ///\code
1699 1707
    ///  b.init();
1700 1708
    ///  for (NodeIt n(gr); n != INVALID; ++n) {
1701 1709
    ///    if (!b.reached(n)) {
1702 1710
    ///      b.addSource(n);
1703 1711
    ///      b.start();
1704 1712
    ///    }
1705 1713
    ///  }
1706 1714
    ///\endcode
1707 1715
    void run() {
1708 1716
      init();
1709 1717
      for (NodeIt it(*_digraph); it != INVALID; ++it) {
1710 1718
        if (!reached(it)) {
1711 1719
          addSource(it);
1712 1720
          start();
1713 1721
        }
1714 1722
      }
1715 1723
    }
1716 1724

	
1717 1725
    ///@}
1718 1726

	
1719 1727
    /// \name Query Functions
1720 1728
    /// The results of the BFS algorithm can be obtained using these
1721 1729
    /// functions.\n
1722 1730
    /// Either \ref run(Node) "run()" or \ref start() should be called
1723 1731
    /// before using them.
1724 1732

	
1725 1733
    ///@{
1726 1734

	
1727 1735
    /// \brief Checks if the given node is reached from the root(s).
1728 1736
    ///
1729 1737
    /// Returns \c true if \c v is reached from the root(s).
1730 1738
    ///
1731 1739
    /// \pre Either \ref run(Node) "run()" or \ref init()
1732 1740
    /// must be called before using this function.
1733 1741
    bool reached(Node v) const { return (*_reached)[v]; }
1734 1742

	
1735 1743
    ///@}
1736 1744

	
1737 1745
  };
1738 1746

	
1739 1747
} //END OF NAMESPACE LEMON
1740 1748

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

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

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

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

	
32 32
namespace lemon {
33 33

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

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

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

	
64 64
  /// \brief Implementation of the Capacity Scaling algorithm for
65 65
  /// finding a \ref min_cost_flow "minimum cost flow".
66 66
  ///
67 67
  /// \ref CapacityScaling implements the capacity scaling version
68 68
  /// of the successive shortest path algorithm for finding a
69 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
  /// and supply values in the algorithm. By default it is \c int.
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
  /// algorithm. By default it is the same as \c V.
82
  /// algorithm. By default, it is the same as \c V.
83
  /// \tparam TR The traits class that defines various types used by the
84
  /// algorithm. By default, it is \ref CapacityScalingDefaultTraits
85
  /// "CapacityScalingDefaultTraits<GR, V, C>".
86
  /// In most cases, this parameter should not be set directly,
87
  /// consider to use the named template parameters instead.
83 88
  ///
84 89
  /// \warning Both number types must be signed and all input data must
85 90
  /// be integer.
86 91
  /// \warning This algorithm does not support negative costs for such
87 92
  /// arcs that have infinite upper bound.
88 93
#ifdef DOXYGEN
89 94
  template <typename GR, typename V, typename C, typename TR>
90 95
#else
91 96
  template < typename GR, typename V = int, typename C = V,
92 97
             typename TR = CapacityScalingDefaultTraits<GR, V, C> >
93 98
#endif
94 99
  class CapacityScaling
95 100
  {
96 101
  public:
97 102

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

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

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

	
111 116
  public:
112 117

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

	
134 139
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
135 140

	
136 141
    typedef std::vector<int> IntVector;
137 142
    typedef std::vector<char> BoolVector;
138 143
    typedef std::vector<Value> ValueVector;
139 144
    typedef std::vector<Cost> CostVector;
140 145

	
141 146
  private:
142 147

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

	
150 155
    // Parameters of the problem
151 156
    bool _have_lower;
152 157
    Value _sum_supply;
153 158

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

	
164 169
    // Node and arc data
165 170
    ValueVector _lower;
166 171
    ValueVector _upper;
167 172
    CostVector _cost;
168 173
    ValueVector _supply;
169 174

	
170 175
    ValueVector _res_cap;
171 176
    CostVector _pi;
172 177
    ValueVector _excess;
173 178
    IntVector _excess_nodes;
174 179
    IntVector _deficit_nodes;
175 180

	
176 181
    Value _delta;
177 182
    int _factor;
178 183
    IntVector _pred;
179 184

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

	
189 194
  private:
190 195

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

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

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

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

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

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

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

	
267 272
        return t;
268 273
      }
269 274

	
270 275
    }; //class ResidualDijkstra
271 276

	
272 277
  public:
273 278

	
274 279
    /// \name Named Template Parameters
275 280
    /// @{
276 281

	
277 282
    template <typename T>
278 283
    struct SetHeapTraits : public Traits {
279 284
      typedef T Heap;
280 285
    };
281 286

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

	
296 301
    /// @}
297 302

	
298 303
  public:
299 304

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

	
317 322
      // Reset data structures
318 323
      reset();
319 324
    }
320 325

	
321 326
    /// \name Parameters
322 327
    /// The parameters of the algorithm can be specified using these
323 328
    /// functions.
324 329

	
325 330
    /// @{
326 331

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

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

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

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

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

	
435 440
    /// \name Execution control
436 441
    /// The algorithm can be executed using \ref run().
437 442

	
438 443
    /// @{
439 444

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

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

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

	
558 563
      _first_out.resize(_node_num + 1);
559 564
      _forward.resize(_res_arc_num);
560 565
      _source.resize(_res_arc_num);
561 566
      _target.resize(_res_arc_num);
562 567
      _reverse.resize(_res_arc_num);
563 568

	
564 569
      _lower.resize(_res_arc_num);
565 570
      _upper.resize(_res_arc_num);
566 571
      _cost.resize(_res_arc_num);
567 572
      _supply.resize(_node_num);
568 573
      
569 574
      _res_cap.resize(_res_arc_num);
570 575
      _pi.resize(_node_num);
571 576
      _excess.resize(_node_num);
572 577
      _pred.resize(_node_num);
573 578

	
574 579
      // Copy the graph
575 580
      int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
576 581
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
577 582
        _node_id[n] = i;
578 583
      }
579 584
      i = 0;
580 585
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
581 586
        _first_out[i] = j;
582 587
        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
583 588
          _arc_idf[a] = j;
584 589
          _forward[j] = true;
585 590
          _source[j] = i;
586 591
          _target[j] = _node_id[_graph.runningNode(a)];
587 592
        }
588 593
        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
589 594
          _arc_idb[a] = j;
590 595
          _forward[j] = false;
591 596
          _source[j] = i;
592 597
          _target[j] = _node_id[_graph.runningNode(a)];
593 598
        }
594 599
        _forward[j] = false;
595 600
        _source[j] = i;
596 601
        _target[j] = _root;
597 602
        _reverse[j] = k;
598 603
        _forward[k] = true;
599 604
        _source[k] = _root;
600 605
        _target[k] = i;
601 606
        _reverse[k] = j;
602 607
        ++j; ++k;
603 608
      }
604 609
      _first_out[i] = j;
605 610
      _first_out[_node_num] = k;
606 611
      for (ArcIt a(_graph); a != INVALID; ++a) {
607 612
        int fi = _arc_idf[a];
608 613
        int bi = _arc_idb[a];
609 614
        _reverse[fi] = bi;
610 615
        _reverse[bi] = fi;
611 616
      }
612 617
      
613 618
      // Reset parameters
614 619
      resetParams();
615 620
      return *this;
616 621
    }
617 622

	
618 623
    /// @}
619 624

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

	
625 630
    /// @{
626 631

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

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

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

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

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

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

	
707 712
    /// @}
708 713

	
709 714
  private:
710 715

	
711 716
    // Initialize the algorithm
712 717
    ProblemType init() {
713 718
      if (_node_num <= 1) return INFEASIBLE;
714 719

	
715 720
      // Check the sum of supply values
716 721
      _sum_supply = 0;
717 722
      for (int i = 0; i != _root; ++i) {
718 723
        _sum_supply += _supply[i];
719 724
      }
720 725
      if (_sum_supply > 0) return INFEASIBLE;
721 726
      
722 727
      // Initialize vectors
723 728
      for (int i = 0; i != _root; ++i) {
724 729
        _pi[i] = 0;
725 730
        _excess[i] = _supply[i];
726 731
      }
727 732

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

	
755 760
      // Handle negative costs
756 761
      for (int i = 0; i != _root; ++i) {
757 762
        last_out = _first_out[i+1] - 1;
758 763
        for (int j = _first_out[i]; j != last_out; ++j) {
759 764
          Value rc = _res_cap[j];
760 765
          if (_cost[j] < 0 && rc > 0) {
761 766
            if (rc >= MAX) return UNBOUNDED;
762 767
            _excess[i] -= rc;
763 768
            _excess[_target[j]] += rc;
764 769
            _res_cap[j] = 0;
765 770
            _res_cap[_reverse[j]] += rc;
766 771
          }
767 772
        }
768 773
      }
769 774
      
770 775
      // Handle GEQ supply type
771 776
      if (_sum_supply < 0) {
772 777
        _pi[_root] = 0;
773 778
        _excess[_root] = -_sum_supply;
774 779
        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
775 780
          int ra = _reverse[a];
776 781
          _res_cap[a] = -_sum_supply + 1;
777 782
          _res_cap[ra] = 0;
778 783
          _cost[a] = 0;
779 784
          _cost[ra] = 0;
780 785
        }
781 786
      } else {
782 787
        _pi[_root] = 0;
783 788
        _excess[_root] = 0;
784 789
        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
785 790
          int ra = _reverse[a];
786 791
          _res_cap[a] = 1;
787 792
          _res_cap[ra] = 0;
788 793
          _cost[a] = 0;
789 794
          _cost[ra] = 0;
790 795
        }
791 796
      }
792 797

	
793 798
      // Initialize delta value
794 799
      if (_factor > 1) {
795 800
        // With scaling
796 801
        Value max_sup = 0, max_dem = 0;
797 802
        for (int i = 0; i != _node_num; ++i) {
798 803
          Value ex = _excess[i];
799 804
          if ( ex > max_sup) max_sup =  ex;
800 805
          if (-ex > max_dem) max_dem = -ex;
801 806
        }
802 807
        Value max_cap = 0;
803 808
        for (int j = 0; j != _res_arc_num; ++j) {
804 809
          if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
805 810
        }
806 811
        max_sup = std::min(std::min(max_sup, max_dem), max_cap);
807 812
        for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;
808 813
      } else {
809 814
        // Without scaling
810 815
        _delta = 1;
811 816
      }
812 817

	
813 818
      return OPTIMAL;
814 819
    }
815 820

	
816 821
    ProblemType start() {
817 822
      // Execute the algorithm
818 823
      ProblemType pt;
819 824
      if (_delta > 1)
820 825
        pt = startWithScaling();
821 826
      else
822 827
        pt = startWithoutScaling();
823 828

	
824 829
      // Handle non-zero lower bounds
825 830
      if (_have_lower) {
826 831
        int limit = _first_out[_root];
827 832
        for (int j = 0; j != limit; ++j) {
828 833
          if (!_forward[j]) _res_cap[j] += _lower[j];
829 834
        }
830 835
      }
831 836

	
832 837
      // Shift potentials if necessary
833 838
      Cost pr = _pi[_root];
834 839
      if (_sum_supply < 0 || pr > 0) {
835 840
        for (int i = 0; i != _node_num; ++i) {
836 841
          _pi[i] -= pr;
837 842
        }        
838 843
      }
839 844
      
840 845
      return pt;
841 846
    }
842 847

	
843 848
    // Execute the capacity scaling algorithm
844 849
    ProblemType startWithScaling() {
845 850
      // Perform capacity scaling phases
846 851
      int s, t;
847 852
      ResidualDijkstra _dijkstra(*this);
848 853
      while (true) {
849 854
        // Saturate all arcs not satisfying the optimality condition
850 855
        int last_out;
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_CIRCULATION_H
20 20
#define LEMON_CIRCULATION_H
21 21

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

	
26 26
///\ingroup max_flow
27 27
///\file
28 28
///\brief Push-relabel algorithm for finding a feasible circulation.
29 29
///
30 30
namespace lemon {
31 31

	
32 32
  /// \brief Default traits class of Circulation class.
33 33
  ///
34 34
  /// Default traits class of Circulation class.
35 35
  ///
36 36
  /// \tparam GR Type of the digraph the algorithm runs on.
37 37
  /// \tparam LM The type of the lower bound map.
38 38
  /// \tparam UM The type of the upper bound (capacity) map.
39 39
  /// \tparam SM The type of the supply map.
40 40
  template <typename GR, typename LM,
41 41
            typename UM, typename SM>
42 42
  struct CirculationDefaultTraits {
43 43

	
44 44
    /// \brief The type of the digraph the algorithm runs on.
45 45
    typedef GR Digraph;
46 46

	
47 47
    /// \brief The type of the lower bound map.
48 48
    ///
49 49
    /// The type of the map that stores the lower bounds on the arcs.
50 50
    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
51 51
    typedef LM LowerMap;
52 52

	
53 53
    /// \brief The type of the upper bound (capacity) map.
54 54
    ///
55 55
    /// The type of the map that stores the upper bounds (capacities)
56 56
    /// on the arcs.
57 57
    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
58 58
    typedef UM UpperMap;
59 59

	
60 60
    /// \brief The type of supply map.
61 61
    ///
62 62
    /// The type of the map that stores the signed supply values of the 
63 63
    /// nodes. 
64 64
    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
65 65
    typedef SM SupplyMap;
66 66

	
67 67
    /// \brief The type of the flow and supply values.
68 68
    typedef typename SupplyMap::Value Value;
69 69

	
70 70
    /// \brief The type of the map that stores the flow values.
71 71
    ///
72 72
    /// The type of the map that stores the flow values.
73 73
    /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
74 74
    /// concept.
75 75
#ifdef DOXYGEN
76 76
    typedef GR::ArcMap<Value> FlowMap;
77 77
#else
78 78
    typedef typename Digraph::template ArcMap<Value> FlowMap;
79 79
#endif
80 80

	
81 81
    /// \brief Instantiates a FlowMap.
82 82
    ///
83 83
    /// This function instantiates a \ref FlowMap.
84 84
    /// \param digraph The digraph for which we would like to define
85 85
    /// the flow map.
86 86
    static FlowMap* createFlowMap(const Digraph& digraph) {
87 87
      return new FlowMap(digraph);
88 88
    }
89 89

	
90 90
    /// \brief The elevator type used by the algorithm.
91 91
    ///
92 92
    /// The elevator type used by the algorithm.
93 93
    ///
94 94
    /// \sa Elevator, LinkedElevator
95 95
#ifdef DOXYGEN
96 96
    typedef lemon::Elevator<GR, GR::Node> Elevator;
97 97
#else
98 98
    typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
99 99
#endif
100 100

	
101 101
    /// \brief Instantiates an Elevator.
102 102
    ///
103 103
    /// This function instantiates an \ref Elevator.
104 104
    /// \param digraph The digraph for which we would like to define
105 105
    /// the elevator.
106 106
    /// \param max_level The maximum level of the elevator.
107 107
    static Elevator* createElevator(const Digraph& digraph, int max_level) {
108 108
      return new Elevator(digraph, max_level);
109 109
    }
110 110

	
111 111
    /// \brief The tolerance used by the algorithm
112 112
    ///
113 113
    /// The tolerance used by the algorithm to handle inexact computation.
114 114
    typedef lemon::Tolerance<Value> Tolerance;
115 115

	
116 116
  };
117 117

	
118 118
  /**
119 119
     \brief Push-relabel algorithm for the network circulation problem.
120 120

	
121 121
     \ingroup max_flow
122 122
     This class implements a push-relabel algorithm for the \e network
123 123
     \e circulation problem.
124 124
     It is to find a feasible circulation when lower and upper bounds
125 125
     are given for the flow values on the arcs and lower bounds are
126 126
     given for the difference between the outgoing and incoming flow
127 127
     at the nodes.
128 128

	
129 129
     The exact formulation of this problem is the following.
130 130
     Let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$
131 131
     \f$upper: A\rightarrow\mathbf{R}\cup\{\infty\}\f$ denote the lower and
132 132
     upper bounds on the arcs, for which \f$lower(uv) \leq upper(uv)\f$
133 133
     holds for all \f$uv\in A\f$, and \f$sup: V\rightarrow\mathbf{R}\f$
134 134
     denotes the signed supply values of the nodes.
135 135
     If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
136 136
     supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
137 137
     \f$-sup(u)\f$ demand.
138 138
     A feasible circulation is an \f$f: A\rightarrow\mathbf{R}\f$
139 139
     solution of the following problem.
140 140

	
141 141
     \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu)
142 142
     \geq sup(u) \quad \forall u\in V, \f]
143 143
     \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A. \f]
144 144
     
145 145
     The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
146 146
     zero or negative in order to have a feasible solution (since the sum
147 147
     of the expressions on the left-hand side of the inequalities is zero).
148 148
     It means that the total demand must be greater or equal to the total
149 149
     supply and all the supplies have to be carried out from the supply nodes,
150 150
     but there could be demands that are not satisfied.
151 151
     If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
152 152
     constraints have to be satisfied with equality, i.e. all demands
153 153
     have to be satisfied and all supplies have to be used.
154 154
     
155 155
     If you need the opposite inequalities in the supply/demand constraints
156 156
     (i.e. the total demand is less than the total supply and all the demands
157 157
     have to be satisfied while there could be supplies that are not used),
158 158
     then you could easily transform the problem to the above form by reversing
159 159
     the direction of the arcs and taking the negative of the supply values
160 160
     (e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
161 161

	
162 162
     This algorithm either calculates a feasible circulation, or provides
163 163
     a \ref barrier() "barrier", which prooves that a feasible soultion
164 164
     cannot exist.
165 165

	
166 166
     Note that this algorithm also provides a feasible solution for the
167 167
     \ref min_cost_flow "minimum cost flow problem".
168 168

	
169 169
     \tparam GR The type of the digraph the algorithm runs on.
170 170
     \tparam LM The type of the lower bound map. The default
171 171
     map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
172 172
     \tparam UM The type of the upper bound (capacity) map.
173 173
     The default map type is \c LM.
174 174
     \tparam SM The type of the supply map. The default map type is
175 175
     \ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>".
176
     \tparam TR The traits class that defines various types used by the
177
     algorithm. By default, it is \ref CirculationDefaultTraits
178
     "CirculationDefaultTraits<GR, LM, UM, SM>".
179
     In most cases, this parameter should not be set directly,
180
     consider to use the named template parameters instead.
176 181
  */
177 182
#ifdef DOXYGEN
178 183
template< typename GR,
179 184
          typename LM,
180 185
          typename UM,
181 186
          typename SM,
182 187
          typename TR >
183 188
#else
184 189
template< typename GR,
185 190
          typename LM = typename GR::template ArcMap<int>,
186 191
          typename UM = LM,
187 192
          typename SM = typename GR::template NodeMap<typename UM::Value>,
188 193
          typename TR = CirculationDefaultTraits<GR, LM, UM, SM> >
189 194
#endif
190 195
  class Circulation {
191 196
  public:
192 197

	
193 198
    ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
194 199
    typedef TR Traits;
195 200
    ///The type of the digraph the algorithm runs on.
196 201
    typedef typename Traits::Digraph Digraph;
197 202
    ///The type of the flow and supply values.
198 203
    typedef typename Traits::Value Value;
199 204

	
200 205
    ///The type of the lower bound map.
201 206
    typedef typename Traits::LowerMap LowerMap;
202 207
    ///The type of the upper bound (capacity) map.
203 208
    typedef typename Traits::UpperMap UpperMap;
204 209
    ///The type of the supply map.
205 210
    typedef typename Traits::SupplyMap SupplyMap;
206 211
    ///The type of the flow map.
207 212
    typedef typename Traits::FlowMap FlowMap;
208 213

	
209 214
    ///The type of the elevator.
210 215
    typedef typename Traits::Elevator Elevator;
211 216
    ///The type of the tolerance.
212 217
    typedef typename Traits::Tolerance Tolerance;
213 218

	
214 219
  private:
215 220

	
216 221
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
217 222

	
218 223
    const Digraph &_g;
219 224
    int _node_num;
220 225

	
221 226
    const LowerMap *_lo;
222 227
    const UpperMap *_up;
223 228
    const SupplyMap *_supply;
224 229

	
225 230
    FlowMap *_flow;
226 231
    bool _local_flow;
227 232

	
228 233
    Elevator* _level;
229 234
    bool _local_level;
230 235

	
231 236
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
232 237
    ExcessMap* _excess;
233 238

	
234 239
    Tolerance _tol;
235 240
    int _el;
236 241

	
237 242
  public:
238 243

	
239 244
    typedef Circulation Create;
240 245

	
241 246
    ///\name Named Template Parameters
242 247

	
243 248
    ///@{
244 249

	
245 250
    template <typename T>
246 251
    struct SetFlowMapTraits : public Traits {
247 252
      typedef T FlowMap;
248 253
      static FlowMap *createFlowMap(const Digraph&) {
249 254
        LEMON_ASSERT(false, "FlowMap is not initialized");
250 255
        return 0; // ignore warnings
251 256
      }
252 257
    };
253 258

	
254 259
    /// \brief \ref named-templ-param "Named parameter" for setting
255 260
    /// FlowMap type
256 261
    ///
257 262
    /// \ref named-templ-param "Named parameter" for setting FlowMap
258 263
    /// type.
259 264
    template <typename T>
260 265
    struct SetFlowMap
261 266
      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
262 267
                           SetFlowMapTraits<T> > {
263 268
      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
264 269
                          SetFlowMapTraits<T> > Create;
265 270
    };
266 271

	
267 272
    template <typename T>
268 273
    struct SetElevatorTraits : public Traits {
269 274
      typedef T Elevator;
270 275
      static Elevator *createElevator(const Digraph&, int) {
271 276
        LEMON_ASSERT(false, "Elevator is not initialized");
272 277
        return 0; // ignore warnings
273 278
      }
274 279
    };
275 280

	
276 281
    /// \brief \ref named-templ-param "Named parameter" for setting
277 282
    /// Elevator type
278 283
    ///
279 284
    /// \ref named-templ-param "Named parameter" for setting Elevator
280 285
    /// type. If this named parameter is used, then an external
281 286
    /// elevator object must be passed to the algorithm using the
282 287
    /// \ref elevator(Elevator&) "elevator()" function before calling
283 288
    /// \ref run() or \ref init().
284 289
    /// \sa SetStandardElevator
285 290
    template <typename T>
286 291
    struct SetElevator
287 292
      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
288 293
                           SetElevatorTraits<T> > {
289 294
      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
290 295
                          SetElevatorTraits<T> > Create;
291 296
    };
292 297

	
293 298
    template <typename T>
294 299
    struct SetStandardElevatorTraits : public Traits {
295 300
      typedef T Elevator;
296 301
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
297 302
        return new Elevator(digraph, max_level);
298 303
      }
299 304
    };
300 305

	
301 306
    /// \brief \ref named-templ-param "Named parameter" for setting
302 307
    /// Elevator type with automatic allocation
303 308
    ///
304 309
    /// \ref named-templ-param "Named parameter" for setting Elevator
305 310
    /// type with automatic allocation.
306 311
    /// The Elevator should have standard constructor interface to be
307 312
    /// able to automatically created by the algorithm (i.e. the
308 313
    /// digraph and the maximum level should be passed to it).
309 314
    /// However, an external elevator object could also be passed to the
310 315
    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
311 316
    /// before calling \ref run() or \ref init().
312 317
    /// \sa SetElevator
313 318
    template <typename T>
314 319
    struct SetStandardElevator
315 320
      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
316 321
                       SetStandardElevatorTraits<T> > {
317 322
      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
318 323
                      SetStandardElevatorTraits<T> > Create;
319 324
    };
320 325

	
321 326
    /// @}
322 327

	
323 328
  protected:
324 329

	
325 330
    Circulation() {}
326 331

	
327 332
  public:
328 333

	
329 334
    /// Constructor.
330 335

	
331 336
    /// The constructor of the class.
332 337
    ///
333 338
    /// \param graph The digraph the algorithm runs on.
334 339
    /// \param lower The lower bounds for the flow values on the arcs.
335 340
    /// \param upper The upper bounds (capacities) for the flow values 
336 341
    /// on the arcs.
337 342
    /// \param supply The signed supply values of the nodes.
338 343
    Circulation(const Digraph &graph, const LowerMap &lower,
339 344
                const UpperMap &upper, const SupplyMap &supply)
340 345
      : _g(graph), _lo(&lower), _up(&upper), _supply(&supply),
341 346
        _flow(NULL), _local_flow(false), _level(NULL), _local_level(false),
342 347
        _excess(NULL) {}
343 348

	
344 349
    /// Destructor.
345 350
    ~Circulation() {
346 351
      destroyStructures();
347 352
    }
348 353

	
349 354

	
350 355
  private:
351 356

	
352 357
    bool checkBoundMaps() {
353 358
      for (ArcIt e(_g);e!=INVALID;++e) {
354 359
        if (_tol.less((*_up)[e], (*_lo)[e])) return false;
355 360
      }
356 361
      return true;
357 362
    }
358 363

	
359 364
    void createStructures() {
360 365
      _node_num = _el = countNodes(_g);
361 366

	
362 367
      if (!_flow) {
363 368
        _flow = Traits::createFlowMap(_g);
364 369
        _local_flow = true;
365 370
      }
366 371
      if (!_level) {
367 372
        _level = Traits::createElevator(_g, _node_num);
368 373
        _local_level = true;
369 374
      }
370 375
      if (!_excess) {
371 376
        _excess = new ExcessMap(_g);
372 377
      }
373 378
    }
374 379

	
375 380
    void destroyStructures() {
376 381
      if (_local_flow) {
377 382
        delete _flow;
378 383
      }
379 384
      if (_local_level) {
380 385
        delete _level;
381 386
      }
382 387
      if (_excess) {
383 388
        delete _excess;
384 389
      }
385 390
    }
386 391

	
387 392
  public:
388 393

	
389 394
    /// Sets the lower bound map.
390 395

	
391 396
    /// Sets the lower bound map.
392 397
    /// \return <tt>(*this)</tt>
393 398
    Circulation& lowerMap(const LowerMap& map) {
394 399
      _lo = &map;
395 400
      return *this;
396 401
    }
397 402

	
398 403
    /// Sets the upper bound (capacity) map.
399 404

	
400 405
    /// Sets the upper bound (capacity) map.
401 406
    /// \return <tt>(*this)</tt>
402 407
    Circulation& upperMap(const UpperMap& map) {
403 408
      _up = &map;
404 409
      return *this;
405 410
    }
406 411

	
407 412
    /// Sets the supply map.
408 413

	
409 414
    /// Sets the supply map.
410 415
    /// \return <tt>(*this)</tt>
411 416
    Circulation& supplyMap(const SupplyMap& map) {
412 417
      _supply = &map;
413 418
      return *this;
414 419
    }
415 420

	
416 421
    /// \brief Sets the flow map.
417 422
    ///
418 423
    /// Sets the flow map.
419 424
    /// If you don't use this function before calling \ref run() or
420 425
    /// \ref init(), an instance will be allocated automatically.
421 426
    /// The destructor deallocates this automatically allocated map,
422 427
    /// of course.
423 428
    /// \return <tt>(*this)</tt>
424 429
    Circulation& flowMap(FlowMap& map) {
425 430
      if (_local_flow) {
426 431
        delete _flow;
427 432
        _local_flow = false;
428 433
      }
429 434
      _flow = &map;
430 435
      return *this;
431 436
    }
432 437

	
433 438
    /// \brief Sets the elevator used by algorithm.
434 439
    ///
435 440
    /// Sets the elevator used by algorithm.
436 441
    /// If you don't use this function before calling \ref run() or
437 442
    /// \ref init(), an instance will be allocated automatically.
438 443
    /// The destructor deallocates this automatically allocated elevator,
439 444
    /// of course.
440 445
    /// \return <tt>(*this)</tt>
441 446
    Circulation& elevator(Elevator& elevator) {
442 447
      if (_local_level) {
443 448
        delete _level;
444 449
        _local_level = false;
445 450
      }
446 451
      _level = &elevator;
447 452
      return *this;
448 453
    }
449 454

	
450 455
    /// \brief Returns a const reference to the elevator.
451 456
    ///
452 457
    /// Returns a const reference to the elevator.
453 458
    ///
454 459
    /// \pre Either \ref run() or \ref init() must be called before
455 460
    /// using this function.
456 461
    const Elevator& elevator() const {
457 462
      return *_level;
458 463
    }
459 464

	
460 465
    /// \brief Sets the tolerance used by the algorithm.
461 466
    ///
462 467
    /// Sets the tolerance object used by the algorithm.
463 468
    /// \return <tt>(*this)</tt>
464 469
    Circulation& tolerance(const Tolerance& tolerance) {
465 470
      _tol = tolerance;
466 471
      return *this;
467 472
    }
468 473

	
469 474
    /// \brief Returns a const reference to the tolerance.
470 475
    ///
471 476
    /// Returns a const reference to the tolerance object used by
472 477
    /// the algorithm.
473 478
    const Tolerance& tolerance() const {
474 479
      return _tol;
475 480
    }
476 481

	
477 482
    /// \name Execution Control
478 483
    /// The simplest way to execute the algorithm is to call \ref run().\n
479 484
    /// If you need better control on the initial solution or the execution,
480 485
    /// you have to call one of the \ref init() functions first, then
481 486
    /// the \ref start() function.
482 487

	
483 488
    ///@{
484 489

	
485 490
    /// Initializes the internal data structures.
486 491

	
487 492
    /// Initializes the internal data structures and sets all flow values
488 493
    /// to the lower bound.
489 494
    void init()
490 495
    {
491 496
      LEMON_DEBUG(checkBoundMaps(),
492 497
        "Upper bounds must be greater or equal to the lower bounds");
493 498

	
494 499
      createStructures();
495 500

	
496 501
      for(NodeIt n(_g);n!=INVALID;++n) {
497 502
        (*_excess)[n] = (*_supply)[n];
498 503
      }
499 504

	
500 505
      for (ArcIt e(_g);e!=INVALID;++e) {
501 506
        _flow->set(e, (*_lo)[e]);
502 507
        (*_excess)[_g.target(e)] += (*_flow)[e];
503 508
        (*_excess)[_g.source(e)] -= (*_flow)[e];
504 509
      }
505 510

	
506 511
      // global relabeling tested, but in general case it provides
507 512
      // worse performance for random digraphs
508 513
      _level->initStart();
509 514
      for(NodeIt n(_g);n!=INVALID;++n)
510 515
        _level->initAddItem(n);
511 516
      _level->initFinish();
512 517
      for(NodeIt n(_g);n!=INVALID;++n)
513 518
        if(_tol.positive((*_excess)[n]))
514 519
          _level->activate(n);
515 520
    }
516 521

	
517 522
    /// Initializes the internal data structures using a greedy approach.
518 523

	
519 524
    /// Initializes the internal data structures using a greedy approach
520 525
    /// to construct the initial solution.
521 526
    void greedyInit()
522 527
    {
523 528
      LEMON_DEBUG(checkBoundMaps(),
524 529
        "Upper bounds must be greater or equal to the lower bounds");
525 530

	
526 531
      createStructures();
527 532

	
528 533
      for(NodeIt n(_g);n!=INVALID;++n) {
529 534
        (*_excess)[n] = (*_supply)[n];
530 535
      }
531 536

	
532 537
      for (ArcIt e(_g);e!=INVALID;++e) {
533 538
        if (!_tol.less(-(*_excess)[_g.target(e)], (*_up)[e])) {
534 539
          _flow->set(e, (*_up)[e]);
535 540
          (*_excess)[_g.target(e)] += (*_up)[e];
536 541
          (*_excess)[_g.source(e)] -= (*_up)[e];
537 542
        } else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) {
538 543
          _flow->set(e, (*_lo)[e]);
539 544
          (*_excess)[_g.target(e)] += (*_lo)[e];
540 545
          (*_excess)[_g.source(e)] -= (*_lo)[e];
541 546
        } else {
542 547
          Value fc = -(*_excess)[_g.target(e)];
543 548
          _flow->set(e, fc);
544 549
          (*_excess)[_g.target(e)] = 0;
545 550
          (*_excess)[_g.source(e)] -= fc;
546 551
        }
547 552
      }
548 553

	
549 554
      _level->initStart();
550 555
      for(NodeIt n(_g);n!=INVALID;++n)
551 556
        _level->initAddItem(n);
552 557
      _level->initFinish();
553 558
      for(NodeIt n(_g);n!=INVALID;++n)
554 559
        if(_tol.positive((*_excess)[n]))
555 560
          _level->activate(n);
556 561
    }
557 562

	
558 563
    ///Executes the algorithm
559 564

	
560 565
    ///This function executes the algorithm.
561 566
    ///
562 567
    ///\return \c true if a feasible circulation is found.
563 568
    ///
564 569
    ///\sa barrier()
565 570
    ///\sa barrierMap()
566 571
    bool start()
567 572
    {
568 573

	
569 574
      Node act;
570 575
      Node bact=INVALID;
571 576
      Node last_activated=INVALID;
572 577
      while((act=_level->highestActive())!=INVALID) {
573 578
        int actlevel=(*_level)[act];
574 579
        int mlevel=_node_num;
575 580
        Value exc=(*_excess)[act];
576 581

	
577 582
        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
578 583
          Node v = _g.target(e);
579 584
          Value fc=(*_up)[e]-(*_flow)[e];
580 585
          if(!_tol.positive(fc)) continue;
581 586
          if((*_level)[v]<actlevel) {
582 587
            if(!_tol.less(fc, exc)) {
583 588
              _flow->set(e, (*_flow)[e] + exc);
584 589
              (*_excess)[v] += exc;
585 590
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
586 591
                _level->activate(v);
587 592
              (*_excess)[act] = 0;
588 593
              _level->deactivate(act);
589 594
              goto next_l;
590 595
            }
591 596
            else {
592 597
              _flow->set(e, (*_up)[e]);
593 598
              (*_excess)[v] += fc;
594 599
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
595 600
                _level->activate(v);
596 601
              exc-=fc;
597 602
            }
598 603
          }
599 604
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
600 605
        }
601 606
        for(InArcIt e(_g,act);e!=INVALID; ++e) {
602 607
          Node v = _g.source(e);
603 608
          Value fc=(*_flow)[e]-(*_lo)[e];
604 609
          if(!_tol.positive(fc)) continue;
605 610
          if((*_level)[v]<actlevel) {
606 611
            if(!_tol.less(fc, exc)) {
607 612
              _flow->set(e, (*_flow)[e] - exc);
608 613
              (*_excess)[v] += exc;
609 614
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
610 615
                _level->activate(v);
611 616
              (*_excess)[act] = 0;
612 617
              _level->deactivate(act);
613 618
              goto next_l;
614 619
            }
615 620
            else {
616 621
              _flow->set(e, (*_lo)[e]);
617 622
              (*_excess)[v] += fc;
618 623
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
619 624
                _level->activate(v);
620 625
              exc-=fc;
621 626
            }
622 627
          }
623 628
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
624 629
        }
625 630

	
626 631
        (*_excess)[act] = exc;
627 632
        if(!_tol.positive(exc)) _level->deactivate(act);
628 633
        else if(mlevel==_node_num) {
629 634
          _level->liftHighestActiveToTop();
630 635
          _el = _node_num;
631 636
          return false;
632 637
        }
633 638
        else {
634 639
          _level->liftHighestActive(mlevel+1);
635 640
          if(_level->onLevel(actlevel)==0) {
636 641
            _el = actlevel;
637 642
            return false;
638 643
          }
639 644
        }
640 645
      next_l:
641 646
        ;
642 647
      }
643 648
      return true;
644 649
    }
645 650

	
646 651
    /// Runs the algorithm.
647 652

	
648 653
    /// This function runs the algorithm.
649 654
    ///
650 655
    /// \return \c true if a feasible circulation is found.
651 656
    ///
652 657
    /// \note Apart from the return value, c.run() is just a shortcut of
653 658
    /// the following code.
654 659
    /// \code
655 660
    ///   c.greedyInit();
656 661
    ///   c.start();
657 662
    /// \endcode
658 663
    bool run() {
659 664
      greedyInit();
660 665
      return start();
661 666
    }
662 667

	
663 668
    /// @}
664 669

	
665 670
    /// \name Query Functions
666 671
    /// The results of the circulation algorithm can be obtained using
667 672
    /// these functions.\n
668 673
    /// Either \ref run() or \ref start() should be called before
669 674
    /// using them.
670 675

	
671 676
    ///@{
672 677

	
673 678
    /// \brief Returns the flow value on the given arc.
674 679
    ///
675 680
    /// Returns the flow value on the given arc.
676 681
    ///
677 682
    /// \pre Either \ref run() or \ref init() must be called before
678 683
    /// using this function.
679 684
    Value flow(const Arc& arc) const {
680 685
      return (*_flow)[arc];
681 686
    }
682 687

	
683 688
    /// \brief Returns a const reference to the flow map.
684 689
    ///
685 690
    /// Returns a const reference to the arc map storing the found flow.
686 691
    ///
687 692
    /// \pre Either \ref run() or \ref init() must be called before
688 693
    /// using this function.
689 694
    const FlowMap& flowMap() const {
690 695
      return *_flow;
691 696
    }
692 697

	
693 698
    /**
694 699
       \brief Returns \c true if the given node is in a barrier.
695 700

	
696 701
       Barrier is a set \e B of nodes for which
697 702

	
698 703
       \f[ \sum_{uv\in A: u\in B} upper(uv) -
699 704
           \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f]
700 705

	
701 706
       holds. The existence of a set with this property prooves that a
702 707
       feasible circualtion cannot exist.
703 708

	
704 709
       This function returns \c true if the given node is in the found
705 710
       barrier. If a feasible circulation is found, the function
706 711
       gives back \c false for every node.
707 712

	
708 713
       \pre Either \ref run() or \ref init() must be called before
709 714
       using this function.
710 715

	
711 716
       \sa barrierMap()
712 717
       \sa checkBarrier()
713 718
    */
714 719
    bool barrier(const Node& node) const
715 720
    {
716 721
      return (*_level)[node] >= _el;
717 722
    }
718 723

	
719 724
    /// \brief Gives back a barrier.
720 725
    ///
721 726
    /// This function sets \c bar to the characteristic vector of the
722 727
    /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
723 728
    /// node map with \c bool (or convertible) value type.
724 729
    ///
725 730
    /// If a feasible circulation is found, the function gives back an
726 731
    /// empty set, so \c bar[v] will be \c false for all nodes \c v.
727 732
    ///
728 733
    /// \note This function calls \ref barrier() for each node,
729 734
    /// so it runs in O(n) time.
730 735
    ///
731 736
    /// \pre Either \ref run() or \ref init() must be called before
732 737
    /// using this function.
733 738
    ///
734 739
    /// \sa barrier()
735 740
    /// \sa checkBarrier()
736 741
    template<class BarrierMap>
737 742
    void barrierMap(BarrierMap &bar) const
738 743
    {
739 744
      for(NodeIt n(_g);n!=INVALID;++n)
740 745
        bar.set(n, (*_level)[n] >= _el);
741 746
    }
742 747

	
743 748
    /// @}
744 749

	
745 750
    /// \name Checker Functions
746 751
    /// The feasibility of the results can be checked using
747 752
    /// these functions.\n
748 753
    /// Either \ref run() or \ref start() should be called before
749 754
    /// using them.
750 755

	
751 756
    ///@{
752 757

	
753 758
    ///Check if the found flow is a feasible circulation
754 759

	
755 760
    ///Check if the found flow is a feasible circulation,
756 761
    ///
757 762
    bool checkFlow() const {
758 763
      for(ArcIt e(_g);e!=INVALID;++e)
759 764
        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
760 765
      for(NodeIt n(_g);n!=INVALID;++n)
761 766
        {
762 767
          Value dif=-(*_supply)[n];
763 768
          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
764 769
          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
765 770
          if(_tol.negative(dif)) return false;
766 771
        }
767 772
      return true;
768 773
    }
769 774

	
770 775
    ///Check whether or not the last execution provides a barrier
771 776

	
772 777
    ///Check whether or not the last execution provides a barrier.
773 778
    ///\sa barrier()
774 779
    ///\sa barrierMap()
775 780
    bool checkBarrier() const
776 781
    {
777 782
      Value delta=0;
778 783
      Value inf_cap = std::numeric_limits<Value>::has_infinity ?
779 784
        std::numeric_limits<Value>::infinity() :
780 785
        std::numeric_limits<Value>::max();
781 786
      for(NodeIt n(_g);n!=INVALID;++n)
782 787
        if(barrier(n))
783 788
          delta-=(*_supply)[n];
784 789
      for(ArcIt e(_g);e!=INVALID;++e)
785 790
        {
786 791
          Node s=_g.source(e);
787 792
          Node t=_g.target(e);
788 793
          if(barrier(s)&&!barrier(t)) {
789 794
            if (_tol.less(inf_cap - (*_up)[e], delta)) return false;
790 795
            delta+=(*_up)[e];
791 796
          }
792 797
          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
793 798
        }
794 799
      return _tol.negative(delta);
795 800
    }
796 801

	
797 802
    /// @}
798 803

	
799 804
  };
800 805

	
801 806
}
802 807

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

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

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

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

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

	
37 37
namespace lemon {
38 38

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

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

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

	
85 85

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

	
89 89
  /// \brief Implementation of the Cost Scaling algorithm for
90 90
  /// finding a \ref min_cost_flow "minimum cost flow".
91 91
  ///
92 92
  /// \ref CostScaling implements a cost scaling algorithm that performs
93 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 100
  /// Most of the parameters of the problem (except for the digraph)
101 101
  /// can be given using separate functions, and the algorithm can be
102 102
  /// executed using the \ref run() function. If some parameters are not
103 103
  /// specified, then default values will be used.
104 104
  ///
105 105
  /// \tparam GR The digraph type the algorithm runs on.
106 106
  /// \tparam V The number type used for flow amounts, capacity bounds
107
  /// and supply values in the algorithm. By default it is \c int.
107
  /// and supply values in the algorithm. By default, it is \c int.
108 108
  /// \tparam C The number type used for costs and potentials in the
109
  /// algorithm. By default it is the same as \c V.
109
  /// algorithm. By default, it is the same as \c V.
110
  /// \tparam TR The traits class that defines various types used by the
111
  /// algorithm. By default, it is \ref CostScalingDefaultTraits
112
  /// "CostScalingDefaultTraits<GR, V, C>".
113
  /// In most cases, this parameter should not be set directly,
114
  /// consider to use the named template parameters instead.
110 115
  ///
111 116
  /// \warning Both number types must be signed and all input data must
112 117
  /// be integer.
113 118
  /// \warning This algorithm does not support negative costs for such
114 119
  /// arcs that have infinite upper bound.
115 120
  ///
116 121
  /// \note %CostScaling provides three different internal methods,
117 122
  /// from which the most efficient one is used by default.
118 123
  /// For more information, see \ref Method.
119 124
#ifdef DOXYGEN
120 125
  template <typename GR, typename V, typename C, typename TR>
121 126
#else
122 127
  template < typename GR, typename V = int, typename C = V,
123 128
             typename TR = CostScalingDefaultTraits<GR, V, C> >
124 129
#endif
125 130
  class CostScaling
126 131
  {
127 132
  public:
128 133

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

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

	
144 148
    /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
145 149
    typedef TR Traits;
146 150

	
147 151
  public:
148 152

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

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

	
195 199
  private:
196 200

	
197 201
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
198 202

	
199 203
    typedef std::vector<int> IntVector;
200 204
    typedef std::vector<char> BoolVector;
201 205
    typedef std::vector<Value> ValueVector;
202 206
    typedef std::vector<Cost> CostVector;
203 207
    typedef std::vector<LargeCost> LargeCostVector;
204 208

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

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

	
227 231
    private:
228 232
      std::vector<Value>& _v;
229 233
    };
230 234

	
231 235
    typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
232 236
    typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
233 237

	
234 238
  private:
235 239

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

	
244 248
    // Parameters of the problem
245 249
    bool _have_lower;
246 250
    Value _sum_supply;
247 251

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

	
258 262
    // Node and arc data
259 263
    ValueVector _lower;
260 264
    ValueVector _upper;
261 265
    CostVector _scost;
262 266
    ValueVector _supply;
263 267

	
264 268
    ValueVector _res_cap;
265 269
    LargeCostVector _cost;
266 270
    LargeCostVector _pi;
267 271
    ValueVector _excess;
268 272
    IntVector _next_out;
269 273
    std::deque<int> _active_nodes;
270 274

	
271 275
    // Data for scaling
272 276
    LargeCost _epsilon;
273 277
    int _alpha;
274 278

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

	
292 296
  public:
293 297

	
294 298
    /// \name Named Template Parameters
295 299
    /// @{
296 300

	
297 301
    template <typename T>
298 302
    struct SetLargeCostTraits : public Traits {
299 303
      typedef T LargeCost;
300 304
    };
301 305

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

	
314 318
    /// @}
315 319

	
316 320
  public:
317 321

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

	
340 344
    /// \name Parameters
341 345
    /// The parameters of the algorithm can be specified using these
342 346
    /// functions.
343 347

	
344 348
    /// @{
345 349

	
346 350
    /// \brief Set the lower bounds on the arcs.
347 351
    ///
348 352
    /// This function sets the lower bounds on the arcs.
349 353
    /// If it is not used before calling \ref run(), the lower bounds
350 354
    /// will be set to zero on all arcs.
351 355
    ///
352 356
    /// \param map An arc map storing the lower bounds.
353 357
    /// Its \c Value type must be convertible to the \c Value type
354 358
    /// of the algorithm.
355 359
    ///
356 360
    /// \return <tt>(*this)</tt>
357 361
    template <typename LowerMap>
358 362
    CostScaling& lowerMap(const LowerMap& map) {
359 363
      _have_lower = true;
360 364
      for (ArcIt a(_graph); a != INVALID; ++a) {
361 365
        _lower[_arc_idf[a]] = map[a];
362 366
        _lower[_arc_idb[a]] = map[a];
363 367
      }
364 368
      return *this;
365 369
    }
366 370

	
367 371
    /// \brief Set the upper bounds (capacities) on the arcs.
368 372
    ///
369 373
    /// This function sets the upper bounds (capacities) on the arcs.
370 374
    /// If it is not used before calling \ref run(), the upper bounds
371 375
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
372 376
    /// unbounded from above).
373 377
    ///
374 378
    /// \param map An arc map storing the upper bounds.
375 379
    /// Its \c Value type must be convertible to the \c Value type
376 380
    /// of the algorithm.
377 381
    ///
378 382
    /// \return <tt>(*this)</tt>
379 383
    template<typename UpperMap>
380 384
    CostScaling& upperMap(const UpperMap& map) {
381 385
      for (ArcIt a(_graph); a != INVALID; ++a) {
382 386
        _upper[_arc_idf[a]] = map[a];
383 387
      }
384 388
      return *this;
385 389
    }
386 390

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

	
407 411
    /// \brief Set the supply values of the nodes.
408 412
    ///
409 413
    /// This function sets the supply values of the nodes.
410 414
    /// If neither this function nor \ref stSupply() is used before
411 415
    /// calling \ref run(), the supply of each node will be set to zero.
412 416
    ///
413 417
    /// \param map A node map storing the supply values.
414 418
    /// Its \c Value type must be convertible to the \c Value type
415 419
    /// of the algorithm.
416 420
    ///
417 421
    /// \return <tt>(*this)</tt>
418 422
    template<typename SupplyMap>
419 423
    CostScaling& supplyMap(const SupplyMap& map) {
420 424
      for (NodeIt n(_graph); n != INVALID; ++n) {
421 425
        _supply[_node_id[n]] = map[n];
422 426
      }
423 427
      return *this;
424 428
    }
425 429

	
426 430
    /// \brief Set single source and target nodes and a supply value.
427 431
    ///
428 432
    /// This function sets a single source node and a single target node
429 433
    /// and the required flow value.
430 434
    /// If neither this function nor \ref supplyMap() is used before
431 435
    /// calling \ref run(), the supply of each node will be set to zero.
432 436
    ///
433 437
    /// Using this function has the same effect as using \ref supplyMap()
434 438
    /// with such a map in which \c k is assigned to \c s, \c -k is
435 439
    /// assigned to \c t and all other nodes have zero supply value.
436 440
    ///
437 441
    /// \param s The source node.
438 442
    /// \param t The target node.
439 443
    /// \param k The required amount of flow from node \c s to node \c t
440 444
    /// (i.e. the supply of \c s and the demand of \c t).
441 445
    ///
442 446
    /// \return <tt>(*this)</tt>
443 447
    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
444 448
      for (int i = 0; i != _res_node_num; ++i) {
445 449
        _supply[i] = 0;
446 450
      }
447 451
      _supply[_node_id[s]] =  k;
448 452
      _supply[_node_id[t]] = -k;
449 453
      return *this;
450 454
    }
451 455
    
452 456
    /// @}
453 457

	
454 458
    /// \name Execution control
455 459
    /// The algorithm can be executed using \ref run().
456 460

	
457 461
    /// @{
458 462

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

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

	
558 562
    /// \brief Reset all the parameters that have been given before.
559 563
    ///
560 564
    /// This function resets all the paramaters that have been given
561 565
    /// before using functions \ref lowerMap(), \ref upperMap(),
562 566
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
563 567
    ///
564 568
    /// It is useful for multiple run() calls. If this function is not
565 569
    /// used, all the parameters given before are kept for the next
566 570
    /// \ref run() call.
567 571
    /// However, the underlying digraph must not be modified after this
568 572
    /// class have been constructed, since it copies and extends the graph.
569 573
    /// \return <tt>(*this)</tt>
570 574
    CostScaling& reset() {
571 575
      // Resize vectors
572 576
      _node_num = countNodes(_graph);
573 577
      _arc_num = countArcs(_graph);
574 578
      _res_node_num = _node_num + 1;
575 579
      _res_arc_num = 2 * (_arc_num + _node_num);
576 580
      _root = _node_num;
577 581

	
578 582
      _first_out.resize(_res_node_num + 1);
579 583
      _forward.resize(_res_arc_num);
580 584
      _source.resize(_res_arc_num);
581 585
      _target.resize(_res_arc_num);
582 586
      _reverse.resize(_res_arc_num);
583 587

	
584 588
      _lower.resize(_res_arc_num);
585 589
      _upper.resize(_res_arc_num);
586 590
      _scost.resize(_res_arc_num);
587 591
      _supply.resize(_res_node_num);
588 592
      
589 593
      _res_cap.resize(_res_arc_num);
590 594
      _cost.resize(_res_arc_num);
591 595
      _pi.resize(_res_node_num);
592 596
      _excess.resize(_res_node_num);
593 597
      _next_out.resize(_res_node_num);
594 598

	
595 599
      _arc_vec.reserve(_res_arc_num);
596 600
      _cost_vec.reserve(_res_arc_num);
597 601

	
598 602
      // Copy the graph
599 603
      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
600 604
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
601 605
        _node_id[n] = i;
602 606
      }
603 607
      i = 0;
604 608
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
605 609
        _first_out[i] = j;
606 610
        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
607 611
          _arc_idf[a] = j;
608 612
          _forward[j] = true;
609 613
          _source[j] = i;
610 614
          _target[j] = _node_id[_graph.runningNode(a)];
611 615
        }
612 616
        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
613 617
          _arc_idb[a] = j;
614 618
          _forward[j] = false;
615 619
          _source[j] = i;
616 620
          _target[j] = _node_id[_graph.runningNode(a)];
617 621
        }
618 622
        _forward[j] = false;
619 623
        _source[j] = i;
620 624
        _target[j] = _root;
621 625
        _reverse[j] = k;
622 626
        _forward[k] = true;
623 627
        _source[k] = _root;
624 628
        _target[k] = i;
625 629
        _reverse[k] = j;
626 630
        ++j; ++k;
627 631
      }
628 632
      _first_out[i] = j;
629 633
      _first_out[_res_node_num] = k;
630 634
      for (ArcIt a(_graph); a != INVALID; ++a) {
631 635
        int fi = _arc_idf[a];
632 636
        int bi = _arc_idb[a];
633 637
        _reverse[fi] = bi;
634 638
        _reverse[bi] = fi;
635 639
      }
636 640
      
637 641
      // Reset parameters
638 642
      resetParams();
639 643
      return *this;
640 644
    }
641 645

	
642 646
    /// @}
643 647

	
644 648
    /// \name Query Functions
645 649
    /// The results of the algorithm can be obtained using these
646 650
    /// functions.\n
647 651
    /// The \ref run() function must be called before using them.
648 652

	
649 653
    /// @{
650 654

	
651 655
    /// \brief Return the total cost of the found flow.
652 656
    ///
653 657
    /// This function returns the total cost of the found flow.
654 658
    /// Its complexity is O(e).
655 659
    ///
656 660
    /// \note The return type of the function can be specified as a
657 661
    /// template parameter. For example,
658 662
    /// \code
659 663
    ///   cs.totalCost<double>();
660 664
    /// \endcode
661 665
    /// It is useful if the total cost cannot be stored in the \c Cost
662 666
    /// type of the algorithm, which is the default return type of the
663 667
    /// function.
664 668
    ///
665 669
    /// \pre \ref run() must be called before using this function.
666 670
    template <typename Number>
667 671
    Number totalCost() const {
668 672
      Number c = 0;
669 673
      for (ArcIt a(_graph); a != INVALID; ++a) {
670 674
        int i = _arc_idb[a];
671 675
        c += static_cast<Number>(_res_cap[i]) *
672 676
             (-static_cast<Number>(_scost[i]));
673 677
      }
674 678
      return c;
675 679
    }
676 680

	
677 681
#ifndef DOXYGEN
678 682
    Cost totalCost() const {
679 683
      return totalCost<Cost>();
680 684
    }
681 685
#endif
682 686

	
683 687
    /// \brief Return the flow on the given arc.
684 688
    ///
685 689
    /// This function returns the flow on the given arc.
686 690
    ///
687 691
    /// \pre \ref run() must be called before using this function.
688 692
    Value flow(const Arc& a) const {
689 693
      return _res_cap[_arc_idb[a]];
690 694
    }
691 695

	
692 696
    /// \brief Return the flow map (the primal solution).
693 697
    ///
694 698
    /// This function copies the flow value on each arc into the given
695 699
    /// map. The \c Value type of the algorithm must be convertible to
696 700
    /// the \c Value type of the map.
697 701
    ///
698 702
    /// \pre \ref run() must be called before using this function.
699 703
    template <typename FlowMap>
700 704
    void flowMap(FlowMap &map) const {
701 705
      for (ArcIt a(_graph); a != INVALID; ++a) {
702 706
        map.set(a, _res_cap[_arc_idb[a]]);
703 707
      }
704 708
    }
705 709

	
706 710
    /// \brief Return the potential (dual value) of the given node.
707 711
    ///
708 712
    /// This function returns the potential (dual value) of the
709 713
    /// given node.
710 714
    ///
711 715
    /// \pre \ref run() must be called before using this function.
712 716
    Cost potential(const Node& n) const {
713 717
      return static_cast<Cost>(_pi[_node_id[n]]);
714 718
    }
715 719

	
716 720
    /// \brief Return the potential map (the dual solution).
717 721
    ///
718 722
    /// This function copies the potential (dual value) of each node
719 723
    /// into the given map.
720 724
    /// The \c Cost type of the algorithm must be convertible to the
721 725
    /// \c Value type of the map.
722 726
    ///
723 727
    /// \pre \ref run() must be called before using this function.
724 728
    template <typename PotentialMap>
725 729
    void potentialMap(PotentialMap &map) const {
726 730
      for (NodeIt n(_graph); n != INVALID; ++n) {
727 731
        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
728 732
      }
729 733
    }
730 734

	
731 735
    /// @}
732 736

	
733 737
  private:
734 738

	
735 739
    // Initialize the algorithm
736 740
    ProblemType init() {
737 741
      if (_res_node_num <= 1) return INFEASIBLE;
738 742

	
739 743
      // Check the sum of supply values
740 744
      _sum_supply = 0;
741 745
      for (int i = 0; i != _root; ++i) {
742 746
        _sum_supply += _supply[i];
743 747
      }
744 748
      if (_sum_supply > 0) return INFEASIBLE;
745 749
      
746 750

	
747 751
      // Initialize vectors
748 752
      for (int i = 0; i != _res_node_num; ++i) {
749 753
        _pi[i] = 0;
750 754
        _excess[i] = _supply[i];
751 755
      }
752 756
      
753 757
      // Remove infinite upper bounds and check negative arcs
754 758
      const Value MAX = std::numeric_limits<Value>::max();
755 759
      int last_out;
756 760
      if (_have_lower) {
757 761
        for (int i = 0; i != _root; ++i) {
758 762
          last_out = _first_out[i+1];
759 763
          for (int j = _first_out[i]; j != last_out; ++j) {
760 764
            if (_forward[j]) {
761 765
              Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
762 766
              if (c >= MAX) return UNBOUNDED;
763 767
              _excess[i] -= c;
764 768
              _excess[_target[j]] += c;
765 769
            }
766 770
          }
767 771
        }
768 772
      } else {
769 773
        for (int i = 0; i != _root; ++i) {
770 774
          last_out = _first_out[i+1];
771 775
          for (int j = _first_out[i]; j != last_out; ++j) {
772 776
            if (_forward[j] && _scost[j] < 0) {
773 777
              Value c = _upper[j];
774 778
              if (c >= MAX) return UNBOUNDED;
775 779
              _excess[i] -= c;
776 780
              _excess[_target[j]] += c;
777 781
            }
778 782
          }
779 783
        }
780 784
      }
781 785
      Value ex, max_cap = 0;
782 786
      for (int i = 0; i != _res_node_num; ++i) {
783 787
        ex = _excess[i];
784 788
        _excess[i] = 0;
785 789
        if (ex < 0) max_cap -= ex;
786 790
      }
787 791
      for (int j = 0; j != _res_arc_num; ++j) {
788 792
        if (_upper[j] >= MAX) _upper[j] = max_cap;
789 793
      }
790 794

	
791 795
      // Initialize the large cost vector and the epsilon parameter
792 796
      _epsilon = 0;
793 797
      LargeCost lc;
794 798
      for (int i = 0; i != _root; ++i) {
795 799
        last_out = _first_out[i+1];
796 800
        for (int j = _first_out[i]; j != last_out; ++j) {
797 801
          lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
798 802
          _cost[j] = lc;
799 803
          if (lc > _epsilon) _epsilon = lc;
800 804
        }
801 805
      }
802 806
      _epsilon /= _alpha;
803 807

	
804 808
      // Initialize maps for Circulation and remove non-zero lower bounds
805 809
      ConstMap<Arc, Value> low(0);
806 810
      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
807 811
      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
808 812
      ValueArcMap cap(_graph), flow(_graph);
809 813
      ValueNodeMap sup(_graph);
810 814
      for (NodeIt n(_graph); n != INVALID; ++n) {
811 815
        sup[n] = _supply[_node_id[n]];
812 816
      }
813 817
      if (_have_lower) {
814 818
        for (ArcIt a(_graph); a != INVALID; ++a) {
815 819
          int j = _arc_idf[a];
816 820
          Value c = _lower[j];
817 821
          cap[a] = _upper[j] - c;
818 822
          sup[_graph.source(a)] -= c;
819 823
          sup[_graph.target(a)] += c;
820 824
        }
821 825
      } else {
822 826
        for (ArcIt a(_graph); a != INVALID; ++a) {
823 827
          cap[a] = _upper[_arc_idf[a]];
824 828
        }
825 829
      }
826 830

	
827 831
      // Find a feasible flow using Circulation
828 832
      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
829 833
        circ(_graph, low, cap, sup);
830 834
      if (!circ.flowMap(flow).run()) return INFEASIBLE;
831 835

	
832 836
      // Set residual capacities and handle GEQ supply type
833 837
      if (_sum_supply < 0) {
834 838
        for (ArcIt a(_graph); a != INVALID; ++a) {
835 839
          Value fa = flow[a];
836 840
          _res_cap[_arc_idf[a]] = cap[a] - fa;
837 841
          _res_cap[_arc_idb[a]] = fa;
838 842
          sup[_graph.source(a)] -= fa;
839 843
          sup[_graph.target(a)] += fa;
840 844
        }
841 845
        for (NodeIt n(_graph); n != INVALID; ++n) {
842 846
          _excess[_node_id[n]] = sup[n];
843 847
        }
844 848
        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
845 849
          int u = _target[a];
846 850
          int ra = _reverse[a];
847 851
          _res_cap[a] = -_sum_supply + 1;
848 852
          _res_cap[ra] = -_excess[u];
849 853
          _cost[a] = 0;
850 854
          _cost[ra] = 0;
851 855
          _excess[u] = 0;
852 856
        }
853 857
      } else {
854 858
        for (ArcIt a(_graph); a != INVALID; ++a) {
855 859
          Value fa = flow[a];
856 860
          _res_cap[_arc_idf[a]] = cap[a] - fa;
857 861
          _res_cap[_arc_idb[a]] = fa;
858 862
        }
859 863
        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
860 864
          int ra = _reverse[a];
861 865
          _res_cap[a] = 1;
862 866
          _res_cap[ra] = 0;
863 867
          _cost[a] = 0;
864 868
          _cost[ra] = 0;
865 869
        }
866 870
      }
867 871
      
868 872
      return OPTIMAL;
869 873
    }
870 874

	
871 875
    // Execute the algorithm and transform the results
872 876
    void start(Method method) {
873 877
      // Maximum path length for partial augment
874 878
      const int MAX_PATH_LENGTH = 4;
875 879
      
876 880
      // Execute the algorithm
877 881
      switch (method) {
878 882
        case PUSH:
879 883
          startPush();
880 884
          break;
881 885
        case AUGMENT:
882 886
          startAugment();
883 887
          break;
884 888
        case PARTIAL_AUGMENT:
885 889
          startAugment(MAX_PATH_LENGTH);
886 890
          break;
887 891
      }
888 892

	
889 893
      // Compute node potentials for the original costs
890 894
      _arc_vec.clear();
891 895
      _cost_vec.clear();
892 896
      for (int j = 0; j != _res_arc_num; ++j) {
893 897
        if (_res_cap[j] > 0) {
894 898
          _arc_vec.push_back(IntPair(_source[j], _target[j]));
895 899
          _cost_vec.push_back(_scost[j]);
896 900
        }
897 901
      }
898 902
      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
899 903

	
900 904
      typename BellmanFord<StaticDigraph, LargeCostArcMap>
901 905
        ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
902 906
      bf.distMap(_pi_map);
903 907
      bf.init(0);
904 908
      bf.start();
905 909

	
906 910
      // Handle non-zero lower bounds
907 911
      if (_have_lower) {
908 912
        int limit = _first_out[_root];

Changeset was too big and was cut off... Show full diff

0 comments (0 inline)