↑ Collapse diff ↑
Ignore white space 6 line context
1 1
CMAKE_MINIMUM_REQUIRED(VERSION 2.6)
2 2

	
3 3
IF(EXISTS ${CMAKE_SOURCE_DIR}/cmake/version.cmake)
4 4
  INCLUDE(${CMAKE_SOURCE_DIR}/cmake/version.cmake)
5 5
ELSE(EXISTS ${CMAKE_SOURCE_DIR}/cmake/version.cmake)
6 6
  SET(PROJECT_NAME "LEMON")
7 7
  SET(PROJECT_VERSION "hg-tip" CACHE STRING "LEMON version string.")
8 8
ENDIF(EXISTS ${CMAKE_SOURCE_DIR}/cmake/version.cmake)
9 9

	
10 10
PROJECT(${PROJECT_NAME})
11 11

	
12 12
SET(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
13 13

	
14 14
INCLUDE(FindDoxygen)
15 15
INCLUDE(FindGhostscript)
16 16
FIND_PACKAGE(GLPK 4.33)
17 17
FIND_PACKAGE(CPLEX)
18 18
FIND_PACKAGE(COIN)
19 19

	
20
ADD_DEFINITIONS(-DHAVE_CONFIG_H)
21

	
22 20
IF(MSVC)
23 21
  SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4250 /wd4355 /wd4800 /wd4996")
24 22
# Suppressed warnings:
25 23
# C4250: 'class1' : inherits 'class2::member' via dominance
26 24
# C4355: 'this' : used in base member initializer list
27 25
# C4800: 'type' : forcing value to bool 'true' or 'false' (performance warning)
28 26
# C4996: 'function': was declared deprecated
29 27
ENDIF(MSVC)
30 28

	
31
ADD_DEFINITIONS(-DHAVE_CONFIG_H)
32

	
33 29
INCLUDE(CheckTypeSize)
34 30
CHECK_TYPE_SIZE("long long" LEMON_LONG_LONG)
35 31

	
36 32
ENABLE_TESTING()
37 33

	
38 34
ADD_SUBDIRECTORY(lemon)
39 35
IF(${CMAKE_SOURCE_DIR} STREQUAL ${PROJECT_SOURCE_DIR})
40 36
  ADD_SUBDIRECTORY(demo)
41 37
  ADD_SUBDIRECTORY(tools)
42 38
  ADD_SUBDIRECTORY(doc)
43 39
  ADD_SUBDIRECTORY(test)
44 40
ENDIF(${CMAKE_SOURCE_DIR} STREQUAL ${PROJECT_SOURCE_DIR})
45 41

	
46 42
IF(${CMAKE_SOURCE_DIR} STREQUAL ${PROJECT_SOURCE_DIR})
47 43
  IF(WIN32)
48 44
    SET(CPACK_PACKAGE_NAME ${PROJECT_NAME})
49 45
    SET(CPACK_PACKAGE_VENDOR "EGRES")
50 46
    SET(CPACK_PACKAGE_DESCRIPTION_SUMMARY
51 47
      "LEMON - Library of Efficient Models and Optimization in Networks")
52 48
    SET(CPACK_RESOURCE_FILE_LICENSE "${PROJECT_SOURCE_DIR}/LICENSE")
53 49

	
54 50
    SET(CPACK_PACKAGE_VERSION ${PROJECT_VERSION})
55 51

	
56 52
    SET(CPACK_PACKAGE_INSTALL_DIRECTORY
57 53
      "${PROJECT_NAME} ${PROJECT_VERSION}")
58 54
    SET(CPACK_PACKAGE_INSTALL_REGISTRY_KEY
59 55
      "${PROJECT_NAME} ${PROJECT_VERSION}")
60 56

	
61 57
    SET(CPACK_COMPONENTS_ALL headers library html_documentation bin)
62 58

	
63 59
    SET(CPACK_COMPONENT_HEADERS_DISPLAY_NAME "C++ headers")
64 60
    SET(CPACK_COMPONENT_LIBRARY_DISPLAY_NAME "Dynamic-link library")
65 61
    SET(CPACK_COMPONENT_BIN_DISPLAY_NAME "Command line utilities")
66 62
    SET(CPACK_COMPONENT_HTML_DOCUMENTATION_DISPLAY_NAME "HTML documentation")
67 63

	
68 64
    SET(CPACK_COMPONENT_HEADERS_DESCRIPTION
69 65
      "C++ header files")
70 66
    SET(CPACK_COMPONENT_LIBRARY_DESCRIPTION
71 67
      "DLL and import library")
72 68
    SET(CPACK_COMPONENT_BIN_DESCRIPTION
73 69
      "Command line utilities")
74 70
    SET(CPACK_COMPONENT_HTML_DOCUMENTATION_DESCRIPTION
75 71
      "Doxygen generated documentation")
76 72

	
77 73
    SET(CPACK_COMPONENT_HEADERS_DEPENDS library)
78 74

	
79 75
    SET(CPACK_COMPONENT_HEADERS_GROUP "Development")
80 76
    SET(CPACK_COMPONENT_LIBRARY_GROUP "Development")
81 77
    SET(CPACK_COMPONENT_HTML_DOCUMENTATION_GROUP "Documentation")
82 78

	
83 79
    SET(CPACK_COMPONENT_GROUP_DEVELOPMENT_DESCRIPTION
84 80
      "Components needed to develop software using LEMON")
85 81
    SET(CPACK_COMPONENT_GROUP_DOCUMENTATION_DESCRIPTION
86 82
      "Documentation of LEMON")
87 83

	
88 84
    SET(CPACK_ALL_INSTALL_TYPES Full Developer)
89 85

	
90 86
    SET(CPACK_COMPONENT_HEADERS_INSTALL_TYPES Developer Full)
91 87
    SET(CPACK_COMPONENT_LIBRARY_INSTALL_TYPES Developer Full)
92 88
    SET(CPACK_COMPONENT_HTML_DOCUMENTATION_INSTALL_TYPES Full)
93 89

	
94 90
    SET(CPACK_GENERATOR "NSIS")
95 91
    SET(CPACK_NSIS_MUI_ICON "${PROJECT_SOURCE_DIR}/cmake/nsis/lemon.ico")
96 92
    SET(CPACK_NSIS_MUI_UNIICON "${PROJECT_SOURCE_DIR}/cmake/nsis/uninstall.ico")
97 93
    #SET(CPACK_PACKAGE_ICON "${PROJECT_SOURCE_DIR}/cmake/nsis\\\\installer.bmp")
98 94
    SET(CPACK_NSIS_INSTALLED_ICON_NAME "bin\\\\lemon.ico")
99 95
    SET(CPACK_NSIS_DISPLAY_NAME "${CPACK_PACKAGE_INSTALL_DIRECTORY} ${PROJECT_NAME}")
100 96
    SET(CPACK_NSIS_HELP_LINK "http:\\\\\\\\lemon.cs.elte.hu")
101 97
    SET(CPACK_NSIS_URL_INFO_ABOUT "http:\\\\\\\\lemon.cs.elte.hu")
102 98
    SET(CPACK_NSIS_CONTACT "lemon-user@lemon.cs.elte.hu")
103 99
    SET(CPACK_NSIS_CREATE_ICONS_EXTRA "
104 100
      CreateShortCut \\\"$SMPROGRAMS\\\\$STARTMENU_FOLDER\\\\Documentation.lnk\\\" \\\"$INSTDIR\\\\share\\\\doc\\\\index.html\\\"
105 101
      ")
106 102
    SET(CPACK_NSIS_DELETE_ICONS_EXTRA "
107 103
      !insertmacro MUI_STARTMENU_GETFOLDER Application $MUI_TEMP
108 104
      Delete \\\"$SMPROGRAMS\\\\$MUI_TEMP\\\\Documentation.lnk\\\"
109 105
      ")
110 106

	
111 107
    INCLUDE(CPack)
112 108
  ENDIF(WIN32)
113 109
ENDIF(${CMAKE_SOURCE_DIR} STREQUAL ${PROJECT_SOURCE_DIR})
Ignore white space 6 line context
1 1
ACLOCAL_AMFLAGS = -I m4
2 2

	
3 3
AM_CXXFLAGS = $(WARNINGCXXFLAGS)
4 4

	
5 5
AM_CPPFLAGS = -I$(top_srcdir) -I$(top_builddir)
6 6
LDADD = $(top_builddir)/lemon/libemon.la
7 7

	
8 8
EXTRA_DIST = \
9 9
	AUTHORS \
10 10
	LICENSE \
11 11
	m4/lx_check_cplex.m4 \
12 12
	m4/lx_check_glpk.m4 \
13 13
	m4/lx_check_soplex.m4 \
14
	m4/lx_check_clp.m4 \
15
	m4/lx_check_cbc.m4 \
14
	m4/lx_check_coin.m4 \
16 15
	CMakeLists.txt \
17 16
	cmake/FindGhostscript.cmake \
17
	cmake/FindCPLEX.cmake \
18 18
	cmake/FindGLPK.cmake \
19
	cmake/FindCOIN.cmake \
19 20
	cmake/version.cmake.in \
20 21
	cmake/version.cmake \
21 22
	cmake/nsis/lemon.ico \
22 23
	cmake/nsis/uninstall.ico
23 24

	
24 25
pkgconfigdir = $(libdir)/pkgconfig
25 26
lemondir = $(pkgincludedir)
26 27
bitsdir = $(lemondir)/bits
27 28
conceptdir = $(lemondir)/concepts
28 29
pkgconfig_DATA =
29 30
lib_LTLIBRARIES =
30 31
lemon_HEADERS =
31 32
bits_HEADERS =
32 33
concept_HEADERS =
33 34
noinst_HEADERS =
34 35
noinst_PROGRAMS =
35 36
bin_PROGRAMS =
36 37
check_PROGRAMS =
37 38
dist_bin_SCRIPTS =
38 39
TESTS =
39 40
XFAIL_TESTS =
40 41

	
41 42
include lemon/Makefile.am
42 43
include test/Makefile.am
43 44
include doc/Makefile.am
44 45
include tools/Makefile.am
45 46

	
46 47
DIST_SUBDIRS = demo
47 48

	
48 49
demo:
49 50
	$(MAKE) $(AM_MAKEFLAGS) -C demo
50 51

	
51 52
MRPROPERFILES = \
52 53
	aclocal.m4 \
53 54
	config.h.in \
54 55
	config.h.in~ \
55 56
	configure \
56 57
	Makefile.in \
57 58
	build-aux/config.guess \
58 59
	build-aux/config.sub \
59 60
	build-aux/depcomp \
60 61
	build-aux/install-sh \
61 62
	build-aux/ltmain.sh \
62 63
	build-aux/missing \
63 64
	doc/doxygen.log
64 65

	
65 66
mrproper:
66 67
	$(MAKE) $(AM_MAKEFLAGS) maintainer-clean
67 68
	-rm -f $(MRPROPERFILES)
68 69

	
69 70
dist-bz2: dist
70 71
	zcat $(PACKAGE)-$(VERSION).tar.gz | \
71 72
	bzip2 --best -c > $(PACKAGE)-$(VERSION).tar.bz2
72 73

	
73 74
distcheck-bz2: distcheck
74 75
	zcat $(PACKAGE)-$(VERSION).tar.gz | \
75 76
	bzip2 --best -c > $(PACKAGE)-$(VERSION).tar.bz2
76 77

	
77 78
.PHONY: demo mrproper dist-bz2 distcheck-bz2
Ignore white space 6 line context
1 1
SET(COIN_ROOT_DIR "" CACHE PATH "COIN root directory")
2 2

	
3 3
FIND_PATH(COIN_INCLUDE_DIR coin/CoinUtilsConfig.h
4
  PATHS ${COIN_ROOT_DIR}/include)
5

	
6
FIND_LIBRARY(COIN_CBC_LIBRARY libCbc
7
  PATHS ${COIN_ROOT_DIR}/lib)
8
FIND_LIBRARY(COIN_CBC_SOLVER_LIBRARY libCbcSolver
9
  PATHS ${COIN_ROOT_DIR}/lib)
10
FIND_LIBRARY(COIN_CGL_LIBRARY libCgl
11
  PATHS ${COIN_ROOT_DIR}/lib)
12
FIND_LIBRARY(COIN_CLP_LIBRARY libClp
13
  PATHS ${COIN_ROOT_DIR}/lib)
14
FIND_LIBRARY(COIN_COIN_UTILS_LIBRARY libCoinUtils
15
  PATHS ${COIN_ROOT_DIR}/lib)
16
FIND_LIBRARY(COIN_OSI_LIBRARY libOsi
17
  PATHS ${COIN_ROOT_DIR}/lib)
18
FIND_LIBRARY(COIN_OSI_CBC_LIBRARY libOsiCbc
19
  PATHS ${COIN_ROOT_DIR}/lib)
20
FIND_LIBRARY(COIN_OSI_CLP_LIBRARY libOsiClp
21
  PATHS ${COIN_ROOT_DIR}/lib)
22
FIND_LIBRARY(COIN_OSI_VOL_LIBRARY libOsiVol
23
  PATHS ${COIN_ROOT_DIR}/lib)
24
FIND_LIBRARY(COIN_VOL_LIBRARY libVol
25
  PATHS ${COIN_ROOT_DIR}/lib)
4
  HINTS ${COIN_ROOT_DIR}/include
5
)
6
FIND_LIBRARY(COIN_CBC_LIBRARY
7
  NAMES Cbc libCbc
8
  HINTS ${COIN_ROOT_DIR}/lib
9
)
10
FIND_LIBRARY(COIN_CBC_SOLVER_LIBRARY
11
  NAMES CbcSolver libCbcSolver
12
  HINTS ${COIN_ROOT_DIR}/lib
13
)
14
FIND_LIBRARY(COIN_CGL_LIBRARY
15
  NAMES Cgl libCgl
16
  HINTS ${COIN_ROOT_DIR}/lib
17
)
18
FIND_LIBRARY(COIN_CLP_LIBRARY
19
  NAMES Clp libClp
20
  HINTS ${COIN_ROOT_DIR}/lib
21
)
22
FIND_LIBRARY(COIN_COIN_UTILS_LIBRARY
23
  NAMES CoinUtils libCoinUtils
24
  HINTS ${COIN_ROOT_DIR}/lib
25
)
26
FIND_LIBRARY(COIN_OSI_LIBRARY
27
  NAMES Osi libOsi
28
  HINTS ${COIN_ROOT_DIR}/lib
29
)
30
FIND_LIBRARY(COIN_OSI_CBC_LIBRARY
31
  NAMES OsiCbc libOsiCbc
32
  HINTS ${COIN_ROOT_DIR}/lib
33
)
34
FIND_LIBRARY(COIN_OSI_CLP_LIBRARY
35
  NAMES OsiClp libOsiClp
36
  HINTS ${COIN_ROOT_DIR}/lib
37
)
38
FIND_LIBRARY(COIN_OSI_VOL_LIBRARY
39
  NAMES OsiVol libOsiVol
40
  HINTS ${COIN_ROOT_DIR}/lib
41
)
42
FIND_LIBRARY(COIN_VOL_LIBRARY
43
  NAMES Vol libVol
44
  HINTS ${COIN_ROOT_DIR}/lib
45
)
26 46

	
27 47
INCLUDE(FindPackageHandleStandardArgs)
28 48
FIND_PACKAGE_HANDLE_STANDARD_ARGS(COIN DEFAULT_MSG
29 49
  COIN_INCLUDE_DIR
30 50
  COIN_CBC_LIBRARY
31 51
  COIN_CBC_SOLVER_LIBRARY
32 52
  COIN_CGL_LIBRARY
33 53
  COIN_CLP_LIBRARY
34 54
  COIN_COIN_UTILS_LIBRARY
35 55
  COIN_OSI_LIBRARY
36 56
  COIN_OSI_CBC_LIBRARY
37 57
  COIN_OSI_CLP_LIBRARY
38 58
  COIN_OSI_VOL_LIBRARY
39 59
  COIN_VOL_LIBRARY
40 60
)
41 61

	
42 62
IF(COIN_FOUND)
43 63
  SET(COIN_INCLUDE_DIRS ${COIN_INCLUDE_DIR})
44 64
  SET(COIN_LIBRARIES "${COIN_CBC_LIBRARY};${COIN_CBC_SOLVER_LIBRARY};${COIN_CGL_LIBRARY};${COIN_CLP_LIBRARY};${COIN_COIN_UTILS_LIBRARY};${COIN_OSI_LIBRARY};${COIN_OSI_CBC_LIBRARY};${COIN_OSI_CLP_LIBRARY};${COIN_OSI_VOL_LIBRARY};${COIN_VOL_LIBRARY}")
45 65
  SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARY};${COIN_COIN_UTILS_LIBRARY}")
46 66
  SET(COIN_CBC_LIBRARIES ${COIN_LIBRARIES})
47 67
ENDIF(COIN_FOUND)
48 68

	
49 69
MARK_AS_ADVANCED(
50 70
  COIN_INCLUDE_DIR
51 71
  COIN_CBC_LIBRARY
52 72
  COIN_CBC_SOLVER_LIBRARY
53 73
  COIN_CGL_LIBRARY
54 74
  COIN_CLP_LIBRARY
55 75
  COIN_COIN_UTILS_LIBRARY
56 76
  COIN_OSI_LIBRARY
57 77
  COIN_OSI_CBC_LIBRARY
58 78
  COIN_OSI_CLP_LIBRARY
59 79
  COIN_OSI_VOL_LIBRARY
60 80
  COIN_VOL_LIBRARY
61 81
)
62 82

	
63 83
IF(COIN_FOUND)
64 84
  SET(LEMON_HAVE_LP TRUE)
65 85
  SET(LEMON_HAVE_MIP TRUE)
66 86
  SET(LEMON_HAVE_CLP TRUE)
67 87
  SET(LEMON_HAVE_CBC TRUE)
68 88
ENDIF(COIN_FOUND)
Ignore white space 6 line context
1
SET(CPLEX_ROOT_DIR "" CACHE PATH "CPLEX root directory")
2

	
1 3
FIND_PATH(CPLEX_INCLUDE_DIR
2 4
  ilcplex/cplex.h
3
  PATHS "C:/ILOG/CPLEX91/include")
4

	
5
  PATHS "C:/ILOG/CPLEX91/include"
6
  PATHS "/opt/ilog/cplex91/include"
7
  HINTS ${CPLEX_ROOT_DIR}/include
8
)
5 9
FIND_LIBRARY(CPLEX_LIBRARY
6
  NAMES cplex91
7
  PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda")
10
  cplex91
11
  PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda"
12
  PATHS "/opt/ilog/cplex91/bin"
13
  HINTS ${CPLEX_ROOT_DIR}/bin
14
)
8 15

	
9 16
INCLUDE(FindPackageHandleStandardArgs)
10 17
FIND_PACKAGE_HANDLE_STANDARD_ARGS(CPLEX DEFAULT_MSG CPLEX_LIBRARY CPLEX_INCLUDE_DIR)
11 18

	
12 19
FIND_PATH(CPLEX_BIN_DIR
13 20
  cplex91.dll
14
  PATHS "C:/ILOG/CPLEX91/bin/x86_win32")
21
  PATHS "C:/ILOG/CPLEX91/bin/x86_win32"
22
)
15 23

	
16 24
IF(CPLEX_FOUND)
17 25
  SET(CPLEX_INCLUDE_DIRS ${CPLEX_INCLUDE_DIR})
18 26
  SET(CPLEX_LIBRARIES ${CPLEX_LIBRARY})
27
  IF(CMAKE_SYSTEM_NAME STREQUAL "Linux")
28
    SET(CPLEX_LIBRARIES "${CPLEX_LIBRARIES};m;pthread")
29
  ENDIF(CMAKE_SYSTEM_NAME STREQUAL "Linux")
19 30
ENDIF(CPLEX_FOUND)
20 31

	
21 32
MARK_AS_ADVANCED(CPLEX_LIBRARY CPLEX_INCLUDE_DIR CPLEX_BIN_DIR)
22 33

	
23 34
IF(CPLEX_FOUND)
24 35
  SET(LEMON_HAVE_LP TRUE)
25 36
  SET(LEMON_HAVE_MIP TRUE)
26 37
  SET(LEMON_HAVE_CPLEX TRUE)
27 38
ENDIF(CPLEX_FOUND)
Ignore white space 6 line context
1
SET(GLPK_ROOT_DIR "" CACHE PATH "GLPK root directory")
2

	
1 3
SET(GLPK_REGKEY "[HKEY_LOCAL_MACHINE\\SOFTWARE\\GnuWin32\\Glpk;InstallPath]")
2 4
GET_FILENAME_COMPONENT(GLPK_ROOT_PATH ${GLPK_REGKEY} ABSOLUTE)
3 5

	
4 6
FIND_PATH(GLPK_INCLUDE_DIR
5 7
  glpk.h
6
  PATHS ${GLPK_REGKEY}/include)
8
  PATHS ${GLPK_REGKEY}/include
9
  HINTS ${GLPK_ROOT_DIR}/include
10
)
11
FIND_LIBRARY(GLPK_LIBRARY
12
  glpk
13
  PATHS ${GLPK_REGKEY}/lib
14
  HINTS ${GLPK_ROOT_DIR}/lib
15
)
7 16

	
8
FIND_LIBRARY(GLPK_LIBRARY
9
  NAMES glpk
10
  PATHS ${GLPK_REGKEY}/lib)
17
IF(GLPK_INCLUDE_DIR AND GLPK_LIBRARY)
18
  FILE(READ ${GLPK_INCLUDE_DIR}/glpk.h GLPK_GLPK_H)
19

	
20
  STRING(REGEX MATCH "define[ ]+GLP_MAJOR_VERSION[ ]+[0-9]+" GLPK_MAJOR_VERSION_LINE "${GLPK_GLPK_H}")
21
  STRING(REGEX REPLACE "define[ ]+GLP_MAJOR_VERSION[ ]+([0-9]+)" "\\1" GLPK_VERSION_MAJOR "${GLPK_MAJOR_VERSION_LINE}")
22

	
23
  STRING(REGEX MATCH "define[ ]+GLP_MINOR_VERSION[ ]+[0-9]+" GLPK_MINOR_VERSION_LINE "${GLPK_GLPK_H}")
24
  STRING(REGEX REPLACE "define[ ]+GLP_MINOR_VERSION[ ]+([0-9]+)" "\\1" GLPK_VERSION_MINOR "${GLPK_MINOR_VERSION_LINE}")
25

	
26
  SET(GLPK_VERSION_STRING "${GLPK_VERSION_MAJOR}.${GLPK_VERSION_MINOR}")
27

	
28
  IF(GLPK_FIND_VERSION)
29
    IF(GLPK_FIND_VERSION_COUNT GREATER 2)
30
      MESSAGE(SEND_ERROR "unexpected version string")
31
    ENDIF(GLPK_FIND_VERSION_COUNT GREATER 2)
32

	
33
    MATH(EXPR GLPK_REQUESTED_VERSION "${GLPK_FIND_VERSION_MAJOR}*100 + ${GLPK_FIND_VERSION_MINOR}")
34
    MATH(EXPR GLPK_FOUND_VERSION "${GLPK_VERSION_MAJOR}*100 + ${GLPK_VERSION_MINOR}")
35

	
36
    IF(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION)
37
      SET(GLPK_PROPER_VERSION_FOUND FALSE)
38
    ELSE(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION)
39
      SET(GLPK_PROPER_VERSION_FOUND TRUE)
40
    ENDIF(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION)
41
  ELSE(GLPK_FIND_VERSION)
42
    SET(GLPK_PROPER_VERSION_FOUND TRUE)
43
  ENDIF(GLPK_FIND_VERSION)
44
ENDIF(GLPK_INCLUDE_DIR AND GLPK_LIBRARY)
11 45

	
12 46
INCLUDE(FindPackageHandleStandardArgs)
13
FIND_PACKAGE_HANDLE_STANDARD_ARGS(GLPK DEFAULT_MSG GLPK_LIBRARY GLPK_INCLUDE_DIR)
47
FIND_PACKAGE_HANDLE_STANDARD_ARGS(GLPK DEFAULT_MSG GLPK_LIBRARY GLPK_INCLUDE_DIR GLPK_PROPER_VERSION_FOUND)
14 48

	
15 49
IF(GLPK_FOUND)
16 50
  SET(GLPK_INCLUDE_DIRS ${GLPK_INCLUDE_DIR})
17 51
  SET(GLPK_LIBRARIES ${GLPK_LIBRARY})
18 52
  SET(GLPK_BIN_DIR ${GLPK_ROOT_PATH}/bin)
19 53
ENDIF(GLPK_FOUND)
20 54

	
21 55
MARK_AS_ADVANCED(GLPK_LIBRARY GLPK_INCLUDE_DIR GLPK_BIN_DIR)
22 56

	
23 57
IF(GLPK_FOUND)
24 58
  SET(LEMON_HAVE_LP TRUE)
25 59
  SET(LEMON_HAVE_MIP TRUE)
26 60
  SET(LEMON_HAVE_GLPK TRUE)
27 61
ENDIF(GLPK_FOUND)
Ignore white space 6 line context
... ...
@@ -163,454 +163,454 @@
163 163
/**
164 164
@defgroup graph_maps Graph Maps
165 165
@ingroup maps
166 166
\brief Special graph-related maps.
167 167

	
168 168
This group contains maps that are specifically designed to assign
169 169
values to the nodes and arcs/edges of graphs.
170 170

	
171 171
If you are looking for the standard graph maps (\c NodeMap, \c ArcMap,
172 172
\c EdgeMap), see the \ref graph_concepts "Graph Structure Concepts".
173 173
*/
174 174

	
175 175
/**
176 176
\defgroup map_adaptors Map Adaptors
177 177
\ingroup maps
178 178
\brief Tools to create new maps from existing ones
179 179

	
180 180
This group contains map adaptors that are used to create "implicit"
181 181
maps from other maps.
182 182

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

	
188 188
The typical usage of this classes is passing implicit maps to
189 189
algorithms.  If a function type algorithm is called then the function
190 190
type map adaptors can be used comfortable. For example let's see the
191 191
usage of map adaptors with the \c graphToEps() function.
192 192
\code
193 193
  Color nodeColor(int deg) {
194 194
    if (deg >= 2) {
195 195
      return Color(0.5, 0.0, 0.5);
196 196
    } else if (deg == 1) {
197 197
      return Color(1.0, 0.5, 1.0);
198 198
    } else {
199 199
      return Color(0.0, 0.0, 0.0);
200 200
    }
201 201
  }
202 202

	
203 203
  Digraph::NodeMap<int> degree_map(graph);
204 204

	
205 205
  graphToEps(graph, "graph.eps")
206 206
    .coords(coords).scaleToA4().undirected()
207 207
    .nodeColors(composeMap(functorToMap(nodeColor), degree_map))
208 208
    .run();
209 209
\endcode
210 210
The \c functorToMap() function makes an \c int to \c Color map from the
211 211
\c nodeColor() function. The \c composeMap() compose the \c degree_map
212 212
and the previously created map. The composed map is a proper function to
213 213
get the color of each node.
214 214

	
215 215
The usage with class type algorithms is little bit harder. In this
216 216
case the function type map adaptors can not be used, because the
217 217
function map adaptors give back temporary objects.
218 218
\code
219 219
  Digraph graph;
220 220

	
221 221
  typedef Digraph::ArcMap<double> DoubleArcMap;
222 222
  DoubleArcMap length(graph);
223 223
  DoubleArcMap speed(graph);
224 224

	
225 225
  typedef DivMap<DoubleArcMap, DoubleArcMap> TimeMap;
226 226
  TimeMap time(length, speed);
227 227

	
228 228
  Dijkstra<Digraph, TimeMap> dijkstra(graph, time);
229 229
  dijkstra.run(source, target);
230 230
\endcode
231 231
We have a length map and a maximum speed map on the arcs of a digraph.
232 232
The minimum time to pass the arc can be calculated as the division of
233 233
the two maps which can be done implicitly with the \c DivMap template
234 234
class. We use the implicit minimum time map as the length map of the
235 235
\c Dijkstra algorithm.
236 236
*/
237 237

	
238 238
/**
239 239
@defgroup matrices Matrices
240 240
@ingroup datas
241 241
\brief Two dimensional data storages implemented in LEMON.
242 242

	
243 243
This group contains two dimensional data storages implemented in LEMON.
244 244
*/
245 245

	
246 246
/**
247 247
@defgroup paths Path Structures
248 248
@ingroup datas
249 249
\brief %Path structures implemented in LEMON.
250 250

	
251 251
This group contains the path structures implemented in LEMON.
252 252

	
253 253
LEMON provides flexible data structures to work with paths.
254 254
All of them have similar interfaces and they can be copied easily with
255 255
assignment operators and copy constructors. This makes it easy and
256 256
efficient to have e.g. the Dijkstra algorithm to store its result in
257 257
any kind of path structure.
258 258

	
259 259
\sa lemon::concepts::Path
260 260
*/
261 261

	
262 262
/**
263 263
@defgroup auxdat Auxiliary Data Structures
264 264
@ingroup datas
265 265
\brief Auxiliary data structures implemented in LEMON.
266 266

	
267 267
This group contains some data structures implemented in LEMON in
268 268
order to make it easier to implement combinatorial algorithms.
269 269
*/
270 270

	
271 271
/**
272 272
@defgroup algs Algorithms
273 273
\brief This group contains the several algorithms
274 274
implemented in LEMON.
275 275

	
276 276
This group contains the several algorithms
277 277
implemented in LEMON.
278 278
*/
279 279

	
280 280
/**
281 281
@defgroup search Graph Search
282 282
@ingroup algs
283 283
\brief Common graph search algorithms.
284 284

	
285 285
This group contains the common graph search algorithms, namely
286 286
\e breadth-first \e search (BFS) and \e depth-first \e search (DFS).
287 287
*/
288 288

	
289 289
/**
290 290
@defgroup shortest_path Shortest Path Algorithms
291 291
@ingroup algs
292 292
\brief Algorithms for finding shortest paths.
293 293

	
294 294
This group contains the algorithms for finding shortest paths in digraphs.
295 295

	
296 296
 - \ref Dijkstra algorithm for finding shortest paths from a source node
297 297
   when all arc lengths are non-negative.
298 298
 - \ref BellmanFord "Bellman-Ford" algorithm for finding shortest paths
299 299
   from a source node when arc lenghts can be either positive or negative,
300 300
   but the digraph should not contain directed cycles with negative total
301 301
   length.
302 302
 - \ref FloydWarshall "Floyd-Warshall" and \ref Johnson "Johnson" algorithms
303 303
   for solving the \e all-pairs \e shortest \e paths \e problem when arc
304 304
   lenghts can be either positive or negative, but the digraph should
305 305
   not contain directed cycles with negative total length.
306 306
 - \ref Suurballe A successive shortest path algorithm for finding
307 307
   arc-disjoint paths between two nodes having minimum total length.
308 308
*/
309 309

	
310 310
/**
311 311
@defgroup max_flow Maximum Flow Algorithms
312 312
@ingroup algs
313 313
\brief Algorithms for finding maximum flows.
314 314

	
315 315
This group contains the algorithms for finding maximum flows and
316 316
feasible circulations.
317 317

	
318 318
The \e maximum \e flow \e problem is to find a flow of maximum value between
319 319
a single source and a single target. Formally, there is a \f$G=(V,A)\f$
320 320
digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and
321 321
\f$s, t \in V\f$ source and target nodes.
322 322
A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the
323 323
following optimization problem.
324 324

	
325 325
\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f]
326 326
\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu)
327 327
    \quad \forall u\in V\setminus\{s,t\} \f]
328 328
\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f]
329 329

	
330 330
LEMON contains several algorithms for solving maximum flow problems:
331 331
- \ref EdmondsKarp Edmonds-Karp algorithm.
332 332
- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm.
333 333
- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees.
334 334
- \ref GoldbergTarjan Preflow push-relabel algorithm with dynamic trees.
335 335

	
336 336
In most cases the \ref Preflow "Preflow" algorithm provides the
337 337
fastest method for computing a maximum flow. All implementations
338 338
provides functions to also query the minimum cut, which is the dual
339 339
problem of the maximum flow.
340 340
*/
341 341

	
342 342
/**
343 343
@defgroup min_cost_flow Minimum Cost Flow Algorithms
344 344
@ingroup algs
345 345

	
346 346
\brief Algorithms for finding minimum cost flows and circulations.
347 347

	
348 348
This group contains the algorithms for finding minimum cost flows and
349 349
circulations.
350 350

	
351 351
The \e minimum \e cost \e flow \e problem is to find a feasible flow of
352 352
minimum total cost from a set of supply nodes to a set of demand nodes
353 353
in a network with capacity constraints (lower and upper bounds)
354 354
and arc costs.
355
Formally, let \f$G=(V,A)\f$ be a digraph,
356
\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
355
Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
356
\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
357 357
upper bounds for the flow values on the arcs, for which
358
\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
359
\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
360
on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
358
\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
359
\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
360
on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
361 361
signed supply values of the nodes.
362 362
If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
363 363
supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
364 364
\f$-sup(u)\f$ demand.
365
A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
365
A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
366 366
of the following optimization problem.
367 367

	
368 368
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
369 369
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
370 370
    sup(u) \quad \forall u\in V \f]
371 371
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
372 372

	
373 373
The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
374 374
zero or negative in order to have a feasible solution (since the sum
375 375
of the expressions on the left-hand side of the inequalities is zero).
376 376
It means that the total demand must be greater or equal to the total
377 377
supply and all the supplies have to be carried out from the supply nodes,
378 378
but there could be demands that are not satisfied.
379 379
If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
380 380
constraints have to be satisfied with equality, i.e. all demands
381 381
have to be satisfied and all supplies have to be used.
382 382

	
383 383
If you need the opposite inequalities in the supply/demand constraints
384 384
(i.e. the total demand is less than the total supply and all the demands
385 385
have to be satisfied while there could be supplies that are not used),
386 386
then you could easily transform the problem to the above form by reversing
387 387
the direction of the arcs and taking the negative of the supply values
388 388
(e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
389 389
However \ref NetworkSimplex algorithm also supports this form directly
390 390
for the sake of convenience.
391 391

	
392 392
A feasible solution for this problem can be found using \ref Circulation.
393 393

	
394 394
Note that the above formulation is actually more general than the usual
395 395
definition of the minimum cost flow problem, in which strict equalities
396 396
are required in the supply/demand contraints, i.e.
397 397

	
398 398
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
399 399
    sup(u) \quad \forall u\in V. \f]
400 400

	
401 401
However if the sum of the supply values is zero, then these two problems
402 402
are equivalent. So if you need the equality form, you have to ensure this
403 403
additional contraint for the algorithms.
404 404

	
405 405
The dual solution of the minimum cost flow problem is represented by node 
406 406
potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
407
An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
407
An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
408 408
is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
409 409
node potentials the following \e complementary \e slackness optimality
410 410
conditions hold.
411 411

	
412 412
 - For all \f$uv\in A\f$ arcs:
413 413
   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
414 414
   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
415 415
   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
416
 - For all \f$u\in V\f$:
416
 - For all \f$u\in V\f$ nodes:
417 417
   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
418 418
     then \f$\pi(u)=0\f$.
419 419
 
420 420
Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
421
\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
421
\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
422 422
\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
423 423

	
424
All algorithms provide dual solution (node potentials) as well
424
All algorithms provide dual solution (node potentials) as well,
425 425
if an optimal flow is found.
426 426

	
427 427
LEMON contains several algorithms for solving minimum cost flow problems.
428 428
 - \ref NetworkSimplex Primal Network Simplex algorithm with various
429 429
   pivot strategies.
430 430
 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
431 431
   cost scaling.
432 432
 - \ref CapacityScaling Successive Shortest %Path algorithm with optional
433 433
   capacity scaling.
434 434
 - \ref CancelAndTighten The Cancel and Tighten algorithm.
435 435
 - \ref CycleCanceling Cycle-Canceling algorithms.
436 436

	
437 437
Most of these implementations support the general inequality form of the
438 438
minimum cost flow problem, but CancelAndTighten and CycleCanceling
439 439
only support the equality form due to the primal method they use.
440 440

	
441 441
In general NetworkSimplex is the most efficient implementation,
442 442
but in special cases other algorithms could be faster.
443 443
For example, if the total supply and/or capacities are rather small,
444 444
CapacityScaling is usually the fastest algorithm (without effective scaling).
445 445
*/
446 446

	
447 447
/**
448 448
@defgroup min_cut Minimum Cut Algorithms
449 449
@ingroup algs
450 450

	
451 451
\brief Algorithms for finding minimum cut in graphs.
452 452

	
453 453
This group contains the algorithms for finding minimum cut in graphs.
454 454

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

	
461 461
\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}}
462 462
    \sum_{uv\in A, u\in X, v\not\in X}cap(uv) \f]
463 463

	
464 464
LEMON contains several algorithms related to minimum cut problems:
465 465

	
466 466
- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut
467 467
  in directed graphs.
468 468
- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for
469 469
  calculating minimum cut in undirected graphs.
470 470
- \ref GomoryHu "Gomory-Hu tree computation" for calculating
471 471
  all-pairs minimum cut in undirected graphs.
472 472

	
473 473
If you want to find minimum cut just between two distinict nodes,
474 474
see the \ref max_flow "maximum flow problem".
475 475
*/
476 476

	
477 477
/**
478 478
@defgroup graph_properties Connectivity and Other Graph Properties
479 479
@ingroup algs
480 480
\brief Algorithms for discovering the graph properties
481 481

	
482 482
This group contains the algorithms for discovering the graph properties
483 483
like connectivity, bipartiteness, euler property, simplicity etc.
484 484

	
485 485
\image html edge_biconnected_components.png
486 486
\image latex edge_biconnected_components.eps "bi-edge-connected components" width=\textwidth
487 487
*/
488 488

	
489 489
/**
490 490
@defgroup planar Planarity Embedding and Drawing
491 491
@ingroup algs
492 492
\brief Algorithms for planarity checking, embedding and drawing
493 493

	
494 494
This group contains the algorithms for planarity checking,
495 495
embedding and drawing.
496 496

	
497 497
\image html planar.png
498 498
\image latex planar.eps "Plane graph" width=\textwidth
499 499
*/
500 500

	
501 501
/**
502 502
@defgroup matching Matching Algorithms
503 503
@ingroup algs
504 504
\brief Algorithms for finding matchings in graphs and bipartite graphs.
505 505

	
506 506
This group contains the algorithms for calculating
507 507
matchings in graphs and bipartite graphs. The general matching problem is
508 508
finding a subset of the edges for which each node has at most one incident
509 509
edge.
510 510

	
511 511
There are several different algorithms for calculate matchings in
512 512
graphs.  The matching problems in bipartite graphs are generally
513 513
easier than in general graphs. The goal of the matching optimization
514 514
can be finding maximum cardinality, maximum weight or minimum cost
515 515
matching. The search can be constrained to find perfect or
516 516
maximum cardinality matching.
517 517

	
518 518
The matching algorithms implemented in LEMON:
519 519
- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm
520 520
  for calculating maximum cardinality matching in bipartite graphs.
521 521
- \ref PrBipartiteMatching Push-relabel algorithm
522 522
  for calculating maximum cardinality matching in bipartite graphs.
523 523
- \ref MaxWeightedBipartiteMatching
524 524
  Successive shortest path algorithm for calculating maximum weighted
525 525
  matching and maximum weighted bipartite matching in bipartite graphs.
526 526
- \ref MinCostMaxBipartiteMatching
527 527
  Successive shortest path algorithm for calculating minimum cost maximum
528 528
  matching in bipartite graphs.
529 529
- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating
530 530
  maximum cardinality matching in general graphs.
531 531
- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating
532 532
  maximum weighted matching in general graphs.
533 533
- \ref MaxWeightedPerfectMatching
534 534
  Edmond's blossom shrinking algorithm for calculating maximum weighted
535 535
  perfect matching in general graphs.
536 536

	
537 537
\image html bipartite_matching.png
538 538
\image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth
539 539
*/
540 540

	
541 541
/**
542 542
@defgroup spantree Minimum Spanning Tree Algorithms
543 543
@ingroup algs
544 544
\brief Algorithms for finding a minimum cost spanning tree in a graph.
545 545

	
546 546
This group contains the algorithms for finding a minimum cost spanning
547 547
tree in a graph.
548 548
*/
549 549

	
550 550
/**
551 551
@defgroup auxalg Auxiliary Algorithms
552 552
@ingroup algs
553 553
\brief Auxiliary algorithms implemented in LEMON.
554 554

	
555 555
This group contains some algorithms implemented in LEMON
556 556
in order to make it easier to implement complex algorithms.
557 557
*/
558 558

	
559 559
/**
560 560
@defgroup approx Approximation Algorithms
561 561
@ingroup algs
562 562
\brief Approximation algorithms.
563 563

	
564 564
This group contains the approximation and heuristic algorithms
565 565
implemented in LEMON.
566 566
*/
567 567

	
568 568
/**
569 569
@defgroup gen_opt_group General Optimization Tools
570 570
\brief This group contains some general optimization frameworks
571 571
implemented in LEMON.
572 572

	
573 573
This group contains some general optimization frameworks
574 574
implemented in LEMON.
575 575
*/
576 576

	
577 577
/**
578 578
@defgroup lp_group Lp and Mip Solvers
579 579
@ingroup gen_opt_group
580 580
\brief Lp and Mip solver interfaces for LEMON.
581 581

	
582 582
This group contains Lp and Mip solver interfaces for LEMON. The
583 583
various LP solvers could be used in the same manner with this
584 584
interface.
585 585
*/
586 586

	
587 587
/**
588 588
@defgroup lp_utils Tools for Lp and Mip Solvers
589 589
@ingroup lp_group
590 590
\brief Helper tools to the Lp and Mip solvers.
591 591

	
592 592
This group adds some helper tools to general optimization framework
593 593
implemented in LEMON.
594 594
*/
595 595

	
596 596
/**
597 597
@defgroup metah Metaheuristics
598 598
@ingroup gen_opt_group
599 599
\brief Metaheuristics for LEMON library.
600 600

	
601 601
This group contains some metaheuristic optimization tools.
602 602
*/
603 603

	
604 604
/**
605 605
@defgroup utils Tools and Utilities
606 606
\brief Tools and utilities for programming in LEMON
607 607

	
608 608
Tools and utilities for programming in LEMON.
609 609
*/
610 610

	
611 611
/**
612 612
@defgroup gutils Basic Graph Utilities
613 613
@ingroup utils
614 614
\brief Simple basic graph utilities.
615 615

	
616 616
This group contains some simple basic graph utilities.
Ignore white space 6 line context
1 1
EXTRA_DIST += \
2 2
	lemon/lemon.pc.in \
3 3
	lemon/CMakeLists.txt
4 4

	
5 5
pkgconfig_DATA += lemon/lemon.pc
6 6

	
7 7
lib_LTLIBRARIES += lemon/libemon.la
8 8

	
9 9
lemon_libemon_la_SOURCES = \
10 10
	lemon/arg_parser.cc \
11 11
	lemon/base.cc \
12 12
	lemon/color.cc \
13 13
	lemon/lp_base.cc \
14 14
	lemon/lp_skeleton.cc \
15 15
	lemon/random.cc \
16 16
	lemon/bits/windows.cc
17 17

	
18

	
18
nodist_lemon_HEADERS = lemon/config.h	
19
	
19 20
lemon_libemon_la_CXXFLAGS = \
20 21
	$(AM_CXXFLAGS) \
21 22
	$(GLPK_CFLAGS) \
22 23
	$(CPLEX_CFLAGS) \
23 24
	$(SOPLEX_CXXFLAGS) \
24 25
	$(CLP_CXXFLAGS) \
25 26
	$(CBC_CXXFLAGS)
26 27

	
27 28
lemon_libemon_la_LDFLAGS = \
28 29
	$(GLPK_LIBS) \
29 30
	$(CPLEX_LIBS) \
30 31
	$(SOPLEX_LIBS) \
31 32
	$(CLP_LIBS) \
32 33
	$(CBC_LIBS)
33 34

	
34 35
if HAVE_GLPK
35 36
lemon_libemon_la_SOURCES += lemon/glpk.cc
36 37
endif
37 38

	
38 39
if HAVE_CPLEX
39 40
lemon_libemon_la_SOURCES += lemon/cplex.cc
40 41
endif
41 42

	
42 43
if HAVE_SOPLEX
43 44
lemon_libemon_la_SOURCES += lemon/soplex.cc
44 45
endif
45 46

	
46 47
if HAVE_CLP
47 48
lemon_libemon_la_SOURCES += lemon/clp.cc
48 49
endif
49 50

	
50 51
if HAVE_CBC
51 52
lemon_libemon_la_SOURCES += lemon/cbc.cc
52 53
endif
53 54

	
54 55
lemon_HEADERS += \
55 56
	lemon/adaptors.h \
56 57
	lemon/arg_parser.h \
57 58
	lemon/assert.h \
58 59
	lemon/bfs.h \
59 60
	lemon/bin_heap.h \
61
	lemon/cbc.h \
60 62
	lemon/circulation.h \
61 63
	lemon/clp.h \
62 64
	lemon/color.h \
63 65
	lemon/concept_check.h \
64
	lemon/config.h \
65 66
	lemon/connectivity.h \
66 67
	lemon/counter.h \
67 68
	lemon/core.h \
68 69
	lemon/cplex.h \
69 70
	lemon/dfs.h \
70 71
	lemon/dijkstra.h \
71 72
	lemon/dim2.h \
72 73
	lemon/dimacs.h \
73 74
	lemon/edge_set.h \
74 75
	lemon/elevator.h \
75 76
	lemon/error.h \
76 77
	lemon/euler.h \
77 78
	lemon/full_graph.h \
78 79
	lemon/glpk.h \
79 80
	lemon/gomory_hu.h \
80 81
	lemon/graph_to_eps.h \
81 82
	lemon/grid_graph.h \
82 83
	lemon/hypercube_graph.h \
83 84
	lemon/kruskal.h \
84 85
	lemon/hao_orlin.h \
85 86
	lemon/lgf_reader.h \
86 87
	lemon/lgf_writer.h \
87 88
	lemon/list_graph.h \
88 89
	lemon/lp.h \
89 90
	lemon/lp_base.h \
90 91
	lemon/lp_skeleton.h \
91 92
	lemon/list_graph.h \
92 93
	lemon/maps.h \
93 94
	lemon/matching.h \
94 95
	lemon/math.h \
95 96
	lemon/min_cost_arborescence.h \
96 97
	lemon/nauty_reader.h \
97 98
	lemon/network_simplex.h \
98 99
	lemon/path.h \
99 100
	lemon/preflow.h \
100 101
	lemon/radix_sort.h \
101 102
	lemon/random.h \
102 103
	lemon/smart_graph.h \
103 104
	lemon/soplex.h \
104 105
	lemon/suurballe.h \
105 106
	lemon/time_measure.h \
106 107
	lemon/tolerance.h \
107 108
	lemon/unionfind.h \
108 109
	lemon/bits/windows.h
109 110

	
110 111
bits_HEADERS += \
111 112
	lemon/bits/alteration_notifier.h \
112 113
	lemon/bits/array_map.h \
113 114
	lemon/bits/base_extender.h \
114 115
	lemon/bits/bezier.h \
115 116
	lemon/bits/default_map.h \
116 117
	lemon/bits/edge_set_extender.h \
117 118
	lemon/bits/enable_if.h \
118 119
	lemon/bits/graph_adaptor_extender.h \
119 120
	lemon/bits/graph_extender.h \
120 121
	lemon/bits/map_extender.h \
121 122
	lemon/bits/path_dump.h \
122 123
	lemon/bits/solver_bits.h \
123 124
	lemon/bits/traits.h \
124 125
	lemon/bits/variant.h \
125 126
	lemon/bits/vector_map.h
126 127

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

	
19 19
#ifndef LEMON_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
    /// \brief The type of the flow values.
68
    typedef typename SupplyMap::Value Flow;
67
    /// \brief The type of the flow and supply values.
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
    typedef typename Digraph::template ArcMap<Flow> FlowMap;
75
    typedef typename Digraph::template ArcMap<Value> FlowMap;
76 76

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

	
86 86
    /// \brief The elevator type used by the algorithm.
87 87
    ///
88 88
    /// The elevator type used by the algorithm.
89 89
    ///
90 90
    /// \sa Elevator
91 91
    /// \sa LinkedElevator
92 92
    typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
93 93

	
94 94
    /// \brief Instantiates an Elevator.
95 95
    ///
96 96
    /// This function instantiates an \ref Elevator.
97 97
    /// \param digraph The digraph for which we would like to define
98 98
    /// the elevator.
99 99
    /// \param max_level The maximum level of the elevator.
100 100
    static Elevator* createElevator(const Digraph& digraph, int max_level) {
101 101
      return new Elevator(digraph, max_level);
102 102
    }
103 103

	
104 104
    /// \brief The tolerance used by the algorithm
105 105
    ///
106 106
    /// The tolerance used by the algorithm to handle inexact computation.
107
    typedef lemon::Tolerance<Flow> Tolerance;
107
    typedef lemon::Tolerance<Value> Tolerance;
108 108

	
109 109
  };
110 110

	
111 111
  /**
112 112
     \brief Push-relabel algorithm for the network circulation problem.
113 113

	
114 114
     \ingroup max_flow
115 115
     This class implements a push-relabel algorithm for the \e network
116 116
     \e circulation problem.
117 117
     It is to find a feasible circulation when lower and upper bounds
118 118
     are given for the flow values on the arcs and lower bounds are
119 119
     given for the difference between the outgoing and incoming flow
120 120
     at the nodes.
121 121

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

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

	
155 155
     This algorithm either calculates a feasible circulation, or provides
156 156
     a \ref barrier() "barrier", which prooves that a feasible soultion
157 157
     cannot exist.
158 158

	
159 159
     Note that this algorithm also provides a feasible solution for the
160 160
     \ref min_cost_flow "minimum cost flow problem".
161 161

	
162 162
     \tparam GR The type of the digraph the algorithm runs on.
163 163
     \tparam LM The type of the lower bound map. The default
164 164
     map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
165 165
     \tparam UM The type of the upper bound (capacity) map.
166 166
     The default map type is \c LM.
167 167
     \tparam SM The type of the supply map. The default map type is
168 168
     \ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>".
169 169
  */
170 170
#ifdef DOXYGEN
171 171
template< typename GR,
172 172
          typename LM,
173 173
          typename UM,
174 174
          typename SM,
175 175
          typename TR >
176 176
#else
177 177
template< typename GR,
178 178
          typename LM = typename GR::template ArcMap<int>,
179 179
          typename UM = LM,
180 180
          typename SM = typename GR::template NodeMap<typename UM::Value>,
181 181
          typename TR = CirculationDefaultTraits<GR, LM, UM, SM> >
182 182
#endif
183 183
  class Circulation {
184 184
  public:
185 185

	
186 186
    ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
187 187
    typedef TR Traits;
188 188
    ///The type of the digraph the algorithm runs on.
189 189
    typedef typename Traits::Digraph Digraph;
190
    ///The type of the flow values.
191
    typedef typename Traits::Flow Flow;
190
    ///The type of the flow and supply values.
191
    typedef typename Traits::Value Value;
192 192

	
193 193
    ///The type of the lower bound map.
194 194
    typedef typename Traits::LowerMap LowerMap;
195 195
    ///The type of the upper bound (capacity) map.
196 196
    typedef typename Traits::UpperMap UpperMap;
197 197
    ///The type of the supply map.
198 198
    typedef typename Traits::SupplyMap SupplyMap;
199 199
    ///The type of the flow map.
200 200
    typedef typename Traits::FlowMap FlowMap;
201 201

	
202 202
    ///The type of the elevator.
203 203
    typedef typename Traits::Elevator Elevator;
204 204
    ///The type of the tolerance.
205 205
    typedef typename Traits::Tolerance Tolerance;
206 206

	
207 207
  private:
208 208

	
209 209
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
210 210

	
211 211
    const Digraph &_g;
212 212
    int _node_num;
213 213

	
214 214
    const LowerMap *_lo;
215 215
    const UpperMap *_up;
216 216
    const SupplyMap *_supply;
217 217

	
218 218
    FlowMap *_flow;
219 219
    bool _local_flow;
220 220

	
221 221
    Elevator* _level;
222 222
    bool _local_level;
223 223

	
224
    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
224
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
225 225
    ExcessMap* _excess;
226 226

	
227 227
    Tolerance _tol;
228 228
    int _el;
229 229

	
230 230
  public:
231 231

	
232 232
    typedef Circulation Create;
233 233

	
234 234
    ///\name Named Template Parameters
235 235

	
236 236
    ///@{
237 237

	
238 238
    template <typename T>
239 239
    struct SetFlowMapTraits : public Traits {
240 240
      typedef T FlowMap;
241 241
      static FlowMap *createFlowMap(const Digraph&) {
242 242
        LEMON_ASSERT(false, "FlowMap is not initialized");
243 243
        return 0; // ignore warnings
244 244
      }
245 245
    };
246 246

	
247 247
    /// \brief \ref named-templ-param "Named parameter" for setting
248 248
    /// FlowMap type
249 249
    ///
250 250
    /// \ref named-templ-param "Named parameter" for setting FlowMap
251 251
    /// type.
252 252
    template <typename T>
253 253
    struct SetFlowMap
254 254
      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
255 255
                           SetFlowMapTraits<T> > {
256 256
      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
257 257
                          SetFlowMapTraits<T> > Create;
258 258
    };
259 259

	
260 260
    template <typename T>
261 261
    struct SetElevatorTraits : public Traits {
262 262
      typedef T Elevator;
263 263
      static Elevator *createElevator(const Digraph&, int) {
264 264
        LEMON_ASSERT(false, "Elevator is not initialized");
265 265
        return 0; // ignore warnings
266 266
      }
267 267
    };
268 268

	
269 269
    /// \brief \ref named-templ-param "Named parameter" for setting
270 270
    /// Elevator type
271 271
    ///
272 272
    /// \ref named-templ-param "Named parameter" for setting Elevator
273 273
    /// type. If this named parameter is used, then an external
274 274
    /// elevator object must be passed to the algorithm using the
275 275
    /// \ref elevator(Elevator&) "elevator()" function before calling
276 276
    /// \ref run() or \ref init().
277 277
    /// \sa SetStandardElevator
278 278
    template <typename T>
279 279
    struct SetElevator
280 280
      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
281 281
                           SetElevatorTraits<T> > {
282 282
      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
283 283
                          SetElevatorTraits<T> > Create;
284 284
    };
285 285

	
286 286
    template <typename T>
287 287
    struct SetStandardElevatorTraits : public Traits {
288 288
      typedef T Elevator;
289 289
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
290 290
        return new Elevator(digraph, max_level);
291 291
      }
292 292
    };
293 293

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

	
314 314
    /// @}
315 315

	
316 316
  protected:
317 317

	
318 318
    Circulation() {}
319 319

	
320 320
  public:
321 321

	
322 322
    /// Constructor.
323 323

	
324 324
    /// The constructor of the class.
325 325
    ///
326 326
    /// \param graph The digraph the algorithm runs on.
327 327
    /// \param lower The lower bounds for the flow values on the arcs.
328 328
    /// \param upper The upper bounds (capacities) for the flow values 
329 329
    /// on the arcs.
330 330
    /// \param supply The signed supply values of the nodes.
331 331
    Circulation(const Digraph &graph, const LowerMap &lower,
332 332
                const UpperMap &upper, const SupplyMap &supply)
333 333
      : _g(graph), _lo(&lower), _up(&upper), _supply(&supply),
334 334
        _flow(NULL), _local_flow(false), _level(NULL), _local_level(false),
335 335
        _excess(NULL) {}
336 336

	
337 337
    /// Destructor.
338 338
    ~Circulation() {
339 339
      destroyStructures();
340 340
    }
341 341

	
342 342

	
343 343
  private:
344 344

	
345 345
    bool checkBoundMaps() {
346 346
      for (ArcIt e(_g);e!=INVALID;++e) {
347 347
        if (_tol.less((*_up)[e], (*_lo)[e])) return false;
348 348
      }
349 349
      return true;
350 350
    }
351 351

	
352 352
    void createStructures() {
353 353
      _node_num = _el = countNodes(_g);
354 354

	
355 355
      if (!_flow) {
356 356
        _flow = Traits::createFlowMap(_g);
357 357
        _local_flow = true;
358 358
      }
359 359
      if (!_level) {
360 360
        _level = Traits::createElevator(_g, _node_num);
361 361
        _local_level = true;
362 362
      }
363 363
      if (!_excess) {
364 364
        _excess = new ExcessMap(_g);
365 365
      }
366 366
    }
367 367

	
368 368
    void destroyStructures() {
369 369
      if (_local_flow) {
370 370
        delete _flow;
371 371
      }
372 372
      if (_local_level) {
373 373
        delete _level;
374 374
      }
375 375
      if (_excess) {
376 376
        delete _excess;
377 377
      }
378 378
    }
379 379

	
380 380
  public:
381 381

	
382 382
    /// Sets the lower bound map.
383 383

	
384 384
    /// Sets the lower bound map.
385 385
    /// \return <tt>(*this)</tt>
386 386
    Circulation& lowerMap(const LowerMap& map) {
387 387
      _lo = &map;
388 388
      return *this;
389 389
    }
390 390

	
391 391
    /// Sets the upper bound (capacity) map.
392 392

	
393 393
    /// Sets the upper bound (capacity) map.
394 394
    /// \return <tt>(*this)</tt>
395 395
    Circulation& upperMap(const UpperMap& map) {
396 396
      _up = &map;
397 397
      return *this;
398 398
    }
399 399

	
400 400
    /// Sets the supply map.
401 401

	
402 402
    /// Sets the supply map.
403 403
    /// \return <tt>(*this)</tt>
404 404
    Circulation& supplyMap(const SupplyMap& map) {
405 405
      _supply = &map;
406 406
      return *this;
407 407
    }
408 408

	
409 409
    /// \brief Sets the flow map.
410 410
    ///
411 411
    /// Sets the flow map.
412 412
    /// If you don't use this function before calling \ref run() or
413 413
    /// \ref init(), an instance will be allocated automatically.
414 414
    /// The destructor deallocates this automatically allocated map,
415 415
    /// of course.
416 416
    /// \return <tt>(*this)</tt>
417 417
    Circulation& flowMap(FlowMap& map) {
418 418
      if (_local_flow) {
419 419
        delete _flow;
420 420
        _local_flow = false;
421 421
      }
422 422
      _flow = &map;
423 423
      return *this;
424 424
    }
425 425

	
426 426
    /// \brief Sets the elevator used by algorithm.
427 427
    ///
428 428
    /// Sets the elevator used by algorithm.
429 429
    /// If you don't use this function before calling \ref run() or
430 430
    /// \ref init(), an instance will be allocated automatically.
431 431
    /// The destructor deallocates this automatically allocated elevator,
432 432
    /// of course.
433 433
    /// \return <tt>(*this)</tt>
434 434
    Circulation& elevator(Elevator& elevator) {
435 435
      if (_local_level) {
436 436
        delete _level;
437 437
        _local_level = false;
438 438
      }
439 439
      _level = &elevator;
440 440
      return *this;
441 441
    }
442 442

	
443 443
    /// \brief Returns a const reference to the elevator.
444 444
    ///
445 445
    /// Returns a const reference to the elevator.
446 446
    ///
447 447
    /// \pre Either \ref run() or \ref init() must be called before
448 448
    /// using this function.
449 449
    const Elevator& elevator() const {
450 450
      return *_level;
451 451
    }
452 452

	
453 453
    /// \brief Sets the tolerance used by algorithm.
454 454
    ///
455 455
    /// Sets the tolerance used by algorithm.
456 456
    Circulation& tolerance(const Tolerance& tolerance) const {
457 457
      _tol = tolerance;
458 458
      return *this;
459 459
    }
460 460

	
461 461
    /// \brief Returns a const reference to the tolerance.
462 462
    ///
463 463
    /// Returns a const reference to the tolerance.
464 464
    const Tolerance& tolerance() const {
465 465
      return tolerance;
466 466
    }
467 467

	
468 468
    /// \name Execution Control
469 469
    /// The simplest way to execute the algorithm is to call \ref run().\n
470 470
    /// If you need more control on the initial solution or the execution,
471 471
    /// first you have to call one of the \ref init() functions, then
472 472
    /// the \ref start() function.
473 473

	
474 474
    ///@{
475 475

	
476 476
    /// Initializes the internal data structures.
477 477

	
478 478
    /// Initializes the internal data structures and sets all flow values
479 479
    /// to the lower bound.
480 480
    void init()
481 481
    {
482 482
      LEMON_DEBUG(checkBoundMaps(),
483 483
        "Upper bounds must be greater or equal to the lower bounds");
484 484

	
485 485
      createStructures();
486 486

	
487 487
      for(NodeIt n(_g);n!=INVALID;++n) {
488 488
        (*_excess)[n] = (*_supply)[n];
489 489
      }
490 490

	
491 491
      for (ArcIt e(_g);e!=INVALID;++e) {
492 492
        _flow->set(e, (*_lo)[e]);
493 493
        (*_excess)[_g.target(e)] += (*_flow)[e];
494 494
        (*_excess)[_g.source(e)] -= (*_flow)[e];
495 495
      }
496 496

	
497 497
      // global relabeling tested, but in general case it provides
498 498
      // worse performance for random digraphs
499 499
      _level->initStart();
500 500
      for(NodeIt n(_g);n!=INVALID;++n)
501 501
        _level->initAddItem(n);
502 502
      _level->initFinish();
503 503
      for(NodeIt n(_g);n!=INVALID;++n)
504 504
        if(_tol.positive((*_excess)[n]))
505 505
          _level->activate(n);
506 506
    }
507 507

	
508 508
    /// Initializes the internal data structures using a greedy approach.
509 509

	
510 510
    /// Initializes the internal data structures using a greedy approach
511 511
    /// to construct the initial solution.
512 512
    void greedyInit()
513 513
    {
514 514
      LEMON_DEBUG(checkBoundMaps(),
515 515
        "Upper bounds must be greater or equal to the lower bounds");
516 516

	
517 517
      createStructures();
518 518

	
519 519
      for(NodeIt n(_g);n!=INVALID;++n) {
520 520
        (*_excess)[n] = (*_supply)[n];
521 521
      }
522 522

	
523 523
      for (ArcIt e(_g);e!=INVALID;++e) {
524 524
        if (!_tol.less(-(*_excess)[_g.target(e)], (*_up)[e])) {
525 525
          _flow->set(e, (*_up)[e]);
526 526
          (*_excess)[_g.target(e)] += (*_up)[e];
527 527
          (*_excess)[_g.source(e)] -= (*_up)[e];
528 528
        } else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) {
529 529
          _flow->set(e, (*_lo)[e]);
530 530
          (*_excess)[_g.target(e)] += (*_lo)[e];
531 531
          (*_excess)[_g.source(e)] -= (*_lo)[e];
532 532
        } else {
533
          Flow fc = -(*_excess)[_g.target(e)];
533
          Value fc = -(*_excess)[_g.target(e)];
534 534
          _flow->set(e, fc);
535 535
          (*_excess)[_g.target(e)] = 0;
536 536
          (*_excess)[_g.source(e)] -= fc;
537 537
        }
538 538
      }
539 539

	
540 540
      _level->initStart();
541 541
      for(NodeIt n(_g);n!=INVALID;++n)
542 542
        _level->initAddItem(n);
543 543
      _level->initFinish();
544 544
      for(NodeIt n(_g);n!=INVALID;++n)
545 545
        if(_tol.positive((*_excess)[n]))
546 546
          _level->activate(n);
547 547
    }
548 548

	
549 549
    ///Executes the algorithm
550 550

	
551 551
    ///This function executes the algorithm.
552 552
    ///
553 553
    ///\return \c true if a feasible circulation is found.
554 554
    ///
555 555
    ///\sa barrier()
556 556
    ///\sa barrierMap()
557 557
    bool start()
558 558
    {
559 559

	
560 560
      Node act;
561 561
      Node bact=INVALID;
562 562
      Node last_activated=INVALID;
563 563
      while((act=_level->highestActive())!=INVALID) {
564 564
        int actlevel=(*_level)[act];
565 565
        int mlevel=_node_num;
566
        Flow exc=(*_excess)[act];
566
        Value exc=(*_excess)[act];
567 567

	
568 568
        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
569 569
          Node v = _g.target(e);
570
          Flow fc=(*_up)[e]-(*_flow)[e];
570
          Value fc=(*_up)[e]-(*_flow)[e];
571 571
          if(!_tol.positive(fc)) continue;
572 572
          if((*_level)[v]<actlevel) {
573 573
            if(!_tol.less(fc, exc)) {
574 574
              _flow->set(e, (*_flow)[e] + exc);
575 575
              (*_excess)[v] += exc;
576 576
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
577 577
                _level->activate(v);
578 578
              (*_excess)[act] = 0;
579 579
              _level->deactivate(act);
580 580
              goto next_l;
581 581
            }
582 582
            else {
583 583
              _flow->set(e, (*_up)[e]);
584 584
              (*_excess)[v] += fc;
585 585
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
586 586
                _level->activate(v);
587 587
              exc-=fc;
588 588
            }
589 589
          }
590 590
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
591 591
        }
592 592
        for(InArcIt e(_g,act);e!=INVALID; ++e) {
593 593
          Node v = _g.source(e);
594
          Flow fc=(*_flow)[e]-(*_lo)[e];
594
          Value fc=(*_flow)[e]-(*_lo)[e];
595 595
          if(!_tol.positive(fc)) continue;
596 596
          if((*_level)[v]<actlevel) {
597 597
            if(!_tol.less(fc, exc)) {
598 598
              _flow->set(e, (*_flow)[e] - exc);
599 599
              (*_excess)[v] += exc;
600 600
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
601 601
                _level->activate(v);
602 602
              (*_excess)[act] = 0;
603 603
              _level->deactivate(act);
604 604
              goto next_l;
605 605
            }
606 606
            else {
607 607
              _flow->set(e, (*_lo)[e]);
608 608
              (*_excess)[v] += fc;
609 609
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
610 610
                _level->activate(v);
611 611
              exc-=fc;
612 612
            }
613 613
          }
614 614
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
615 615
        }
616 616

	
617 617
        (*_excess)[act] = exc;
618 618
        if(!_tol.positive(exc)) _level->deactivate(act);
619 619
        else if(mlevel==_node_num) {
620 620
          _level->liftHighestActiveToTop();
621 621
          _el = _node_num;
622 622
          return false;
623 623
        }
624 624
        else {
625 625
          _level->liftHighestActive(mlevel+1);
626 626
          if(_level->onLevel(actlevel)==0) {
627 627
            _el = actlevel;
628 628
            return false;
629 629
          }
630 630
        }
631 631
      next_l:
632 632
        ;
633 633
      }
634 634
      return true;
635 635
    }
636 636

	
637 637
    /// Runs the algorithm.
638 638

	
639 639
    /// This function runs the algorithm.
640 640
    ///
641 641
    /// \return \c true if a feasible circulation is found.
642 642
    ///
643 643
    /// \note Apart from the return value, c.run() is just a shortcut of
644 644
    /// the following code.
645 645
    /// \code
646 646
    ///   c.greedyInit();
647 647
    ///   c.start();
648 648
    /// \endcode
649 649
    bool run() {
650 650
      greedyInit();
651 651
      return start();
652 652
    }
653 653

	
654 654
    /// @}
655 655

	
656 656
    /// \name Query Functions
657 657
    /// The results of the circulation algorithm can be obtained using
658 658
    /// these functions.\n
659 659
    /// Either \ref run() or \ref start() should be called before
660 660
    /// using them.
661 661

	
662 662
    ///@{
663 663

	
664
    /// \brief Returns the flow on the given arc.
664
    /// \brief Returns the flow value on the given arc.
665 665
    ///
666
    /// Returns the flow on the given arc.
666
    /// Returns the flow value on the given arc.
667 667
    ///
668 668
    /// \pre Either \ref run() or \ref init() must be called before
669 669
    /// using this function.
670
    Flow flow(const Arc& arc) const {
670
    Value flow(const Arc& arc) const {
671 671
      return (*_flow)[arc];
672 672
    }
673 673

	
674 674
    /// \brief Returns a const reference to the flow map.
675 675
    ///
676 676
    /// Returns a const reference to the arc map storing the found flow.
677 677
    ///
678 678
    /// \pre Either \ref run() or \ref init() must be called before
679 679
    /// using this function.
680 680
    const FlowMap& flowMap() const {
681 681
      return *_flow;
682 682
    }
683 683

	
684 684
    /**
685 685
       \brief Returns \c true if the given node is in a barrier.
686 686

	
687 687
       Barrier is a set \e B of nodes for which
688 688

	
689 689
       \f[ \sum_{uv\in A: u\in B} upper(uv) -
690 690
           \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f]
691 691

	
692 692
       holds. The existence of a set with this property prooves that a
693 693
       feasible circualtion cannot exist.
694 694

	
695 695
       This function returns \c true if the given node is in the found
696 696
       barrier. If a feasible circulation is found, the function
697 697
       gives back \c false for every node.
698 698

	
699 699
       \pre Either \ref run() or \ref init() must be called before
700 700
       using this function.
701 701

	
702 702
       \sa barrierMap()
703 703
       \sa checkBarrier()
704 704
    */
705 705
    bool barrier(const Node& node) const
706 706
    {
707 707
      return (*_level)[node] >= _el;
708 708
    }
709 709

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

	
734 734
    /// @}
735 735

	
736 736
    /// \name Checker Functions
737 737
    /// The feasibility of the results can be checked using
738 738
    /// these functions.\n
739 739
    /// Either \ref run() or \ref start() should be called before
740 740
    /// using them.
741 741

	
742 742
    ///@{
743 743

	
744 744
    ///Check if the found flow is a feasible circulation
745 745

	
746 746
    ///Check if the found flow is a feasible circulation,
747 747
    ///
748 748
    bool checkFlow() const {
749 749
      for(ArcIt e(_g);e!=INVALID;++e)
750 750
        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
751 751
      for(NodeIt n(_g);n!=INVALID;++n)
752 752
        {
753
          Flow dif=-(*_supply)[n];
753
          Value dif=-(*_supply)[n];
754 754
          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
755 755
          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
756 756
          if(_tol.negative(dif)) return false;
757 757
        }
758 758
      return true;
759 759
    }
760 760

	
761 761
    ///Check whether or not the last execution provides a barrier
762 762

	
763 763
    ///Check whether or not the last execution provides a barrier.
764 764
    ///\sa barrier()
765 765
    ///\sa barrierMap()
766 766
    bool checkBarrier() const
767 767
    {
768
      Flow delta=0;
769
      Flow inf_cap = std::numeric_limits<Flow>::has_infinity ?
770
        std::numeric_limits<Flow>::infinity() :
771
        std::numeric_limits<Flow>::max();
768
      Value delta=0;
769
      Value inf_cap = std::numeric_limits<Value>::has_infinity ?
770
        std::numeric_limits<Value>::infinity() :
771
        std::numeric_limits<Value>::max();
772 772
      for(NodeIt n(_g);n!=INVALID;++n)
773 773
        if(barrier(n))
774 774
          delta-=(*_supply)[n];
775 775
      for(ArcIt e(_g);e!=INVALID;++e)
776 776
        {
777 777
          Node s=_g.source(e);
778 778
          Node t=_g.target(e);
779 779
          if(barrier(s)&&!barrier(t)) {
780 780
            if (_tol.less(inf_cap - (*_up)[e], delta)) return false;
781 781
            delta+=(*_up)[e];
782 782
          }
783 783
          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
784 784
        }
785 785
      return _tol.negative(delta);
786 786
    }
787 787

	
788 788
    /// @}
789 789

	
790 790
  };
791 791

	
792 792
}
793 793

	
794 794
#endif
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_CORE_H
20 20
#define LEMON_CORE_H
21 21

	
22 22
#include <vector>
23 23
#include <algorithm>
24 24

	
25
#include <lemon/core.h>
25
#include <lemon/config.h>
26 26
#include <lemon/bits/enable_if.h>
27 27
#include <lemon/bits/traits.h>
28 28
#include <lemon/assert.h>
29 29

	
30 30
///\file
31 31
///\brief LEMON core utilities.
32 32
///
33 33
///This header file contains core utilities for LEMON.
34 34
///It is automatically included by all graph types, therefore it usually
35 35
///do not have to be included directly.
36 36

	
37 37
namespace lemon {
38 38

	
39 39
  /// \brief Dummy type to make it easier to create invalid iterators.
40 40
  ///
41 41
  /// Dummy type to make it easier to create invalid iterators.
42 42
  /// See \ref INVALID for the usage.
43 43
  struct Invalid {
44 44
  public:
45 45
    bool operator==(Invalid) { return true;  }
46 46
    bool operator!=(Invalid) { return false; }
47 47
    bool operator< (Invalid) { return false; }
48 48
  };
49 49

	
50 50
  /// \brief Invalid iterators.
51 51
  ///
52 52
  /// \ref Invalid is a global type that converts to each iterator
53 53
  /// in such a way that the value of the target iterator will be invalid.
54 54
#ifdef LEMON_ONLY_TEMPLATES
55 55
  const Invalid INVALID = Invalid();
56 56
#else
57 57
  extern const Invalid INVALID;
58 58
#endif
59 59

	
60 60
  /// \addtogroup gutils
61 61
  /// @{
62 62

	
63 63
  ///Create convenience typedefs for the digraph types and iterators
64 64

	
65 65
  ///This \c \#define creates convenient type definitions for the following
66 66
  ///types of \c Digraph: \c Node,  \c NodeIt, \c Arc, \c ArcIt, \c InArcIt,
67 67
  ///\c OutArcIt, \c BoolNodeMap, \c IntNodeMap, \c DoubleNodeMap,
68 68
  ///\c BoolArcMap, \c IntArcMap, \c DoubleArcMap.
69 69
  ///
70 70
  ///\note If the graph type is a dependent type, ie. the graph type depend
71 71
  ///on a template parameter, then use \c TEMPLATE_DIGRAPH_TYPEDEFS()
72 72
  ///macro.
73 73
#define DIGRAPH_TYPEDEFS(Digraph)                                       \
74 74
  typedef Digraph::Node Node;                                           \
75 75
  typedef Digraph::NodeIt NodeIt;                                       \
76 76
  typedef Digraph::Arc Arc;                                             \
77 77
  typedef Digraph::ArcIt ArcIt;                                         \
78 78
  typedef Digraph::InArcIt InArcIt;                                     \
79 79
  typedef Digraph::OutArcIt OutArcIt;                                   \
80 80
  typedef Digraph::NodeMap<bool> BoolNodeMap;                           \
81 81
  typedef Digraph::NodeMap<int> IntNodeMap;                             \
82 82
  typedef Digraph::NodeMap<double> DoubleNodeMap;                       \
83 83
  typedef Digraph::ArcMap<bool> BoolArcMap;                             \
84 84
  typedef Digraph::ArcMap<int> IntArcMap;                               \
85 85
  typedef Digraph::ArcMap<double> DoubleArcMap
86 86

	
87 87
  ///Create convenience typedefs for the digraph types and iterators
88 88

	
89 89
  ///\see DIGRAPH_TYPEDEFS
90 90
  ///
91 91
  ///\note Use this macro, if the graph type is a dependent type,
92 92
  ///ie. the graph type depend on a template parameter.
93 93
#define TEMPLATE_DIGRAPH_TYPEDEFS(Digraph)                              \
94 94
  typedef typename Digraph::Node Node;                                  \
95 95
  typedef typename Digraph::NodeIt NodeIt;                              \
96 96
  typedef typename Digraph::Arc Arc;                                    \
97 97
  typedef typename Digraph::ArcIt ArcIt;                                \
98 98
  typedef typename Digraph::InArcIt InArcIt;                            \
99 99
  typedef typename Digraph::OutArcIt OutArcIt;                          \
100 100
  typedef typename Digraph::template NodeMap<bool> BoolNodeMap;         \
101 101
  typedef typename Digraph::template NodeMap<int> IntNodeMap;           \
102 102
  typedef typename Digraph::template NodeMap<double> DoubleNodeMap;     \
103 103
  typedef typename Digraph::template ArcMap<bool> BoolArcMap;           \
104 104
  typedef typename Digraph::template ArcMap<int> IntArcMap;             \
105 105
  typedef typename Digraph::template ArcMap<double> DoubleArcMap
106 106

	
107 107
  ///Create convenience typedefs for the graph types and iterators
108 108

	
109 109
  ///This \c \#define creates the same convenient type definitions as defined
110 110
  ///by \ref DIGRAPH_TYPEDEFS(Graph) and six more, namely it creates
111 111
  ///\c Edge, \c EdgeIt, \c IncEdgeIt, \c BoolEdgeMap, \c IntEdgeMap,
112 112
  ///\c DoubleEdgeMap.
113 113
  ///
114 114
  ///\note If the graph type is a dependent type, ie. the graph type depend
115 115
  ///on a template parameter, then use \c TEMPLATE_GRAPH_TYPEDEFS()
116 116
  ///macro.
117 117
#define GRAPH_TYPEDEFS(Graph)                                           \
118 118
  DIGRAPH_TYPEDEFS(Graph);                                              \
119 119
  typedef Graph::Edge Edge;                                             \
120 120
  typedef Graph::EdgeIt EdgeIt;                                         \
121 121
  typedef Graph::IncEdgeIt IncEdgeIt;                                   \
122 122
  typedef Graph::EdgeMap<bool> BoolEdgeMap;                             \
123 123
  typedef Graph::EdgeMap<int> IntEdgeMap;                               \
124 124
  typedef Graph::EdgeMap<double> DoubleEdgeMap
125 125

	
126 126
  ///Create convenience typedefs for the graph types and iterators
127 127

	
128 128
  ///\see GRAPH_TYPEDEFS
129 129
  ///
130 130
  ///\note Use this macro, if the graph type is a dependent type,
131 131
  ///ie. the graph type depend on a template parameter.
132 132
#define TEMPLATE_GRAPH_TYPEDEFS(Graph)                                  \
133 133
  TEMPLATE_DIGRAPH_TYPEDEFS(Graph);                                     \
134 134
  typedef typename Graph::Edge Edge;                                    \
135 135
  typedef typename Graph::EdgeIt EdgeIt;                                \
136 136
  typedef typename Graph::IncEdgeIt IncEdgeIt;                          \
137 137
  typedef typename Graph::template EdgeMap<bool> BoolEdgeMap;           \
138 138
  typedef typename Graph::template EdgeMap<int> IntEdgeMap;             \
139 139
  typedef typename Graph::template EdgeMap<double> DoubleEdgeMap
140 140

	
141 141
  /// \brief Function to count the items in a graph.
142 142
  ///
143 143
  /// This function counts the items (nodes, arcs etc.) in a graph.
144 144
  /// The complexity of the function is linear because
145 145
  /// it iterates on all of the items.
146 146
  template <typename Graph, typename Item>
147 147
  inline int countItems(const Graph& g) {
148 148
    typedef typename ItemSetTraits<Graph, Item>::ItemIt ItemIt;
149 149
    int num = 0;
150 150
    for (ItemIt it(g); it != INVALID; ++it) {
151 151
      ++num;
152 152
    }
153 153
    return num;
154 154
  }
155 155

	
156 156
  // Node counting:
157 157

	
158 158
  namespace _core_bits {
159 159

	
160 160
    template <typename Graph, typename Enable = void>
161 161
    struct CountNodesSelector {
162 162
      static int count(const Graph &g) {
163 163
        return countItems<Graph, typename Graph::Node>(g);
164 164
      }
165 165
    };
166 166

	
167 167
    template <typename Graph>
168 168
    struct CountNodesSelector<
169 169
      Graph, typename
170 170
      enable_if<typename Graph::NodeNumTag, void>::type>
171 171
    {
172 172
      static int count(const Graph &g) {
173 173
        return g.nodeNum();
174 174
      }
175 175
    };
176 176
  }
177 177

	
178 178
  /// \brief Function to count the nodes in the graph.
179 179
  ///
180 180
  /// This function counts the nodes in the graph.
181 181
  /// The complexity of the function is <em>O</em>(<em>n</em>), but for some
182 182
  /// graph structures it is specialized to run in <em>O</em>(1).
183 183
  ///
184 184
  /// \note If the graph contains a \c nodeNum() member function and a
185 185
  /// \c NodeNumTag tag then this function calls directly the member
186 186
  /// function to query the cardinality of the node set.
187 187
  template <typename Graph>
188 188
  inline int countNodes(const Graph& g) {
189 189
    return _core_bits::CountNodesSelector<Graph>::count(g);
190 190
  }
191 191

	
192 192
  // Arc counting:
193 193

	
194 194
  namespace _core_bits {
195 195

	
196 196
    template <typename Graph, typename Enable = void>
197 197
    struct CountArcsSelector {
198 198
      static int count(const Graph &g) {
199 199
        return countItems<Graph, typename Graph::Arc>(g);
200 200
      }
201 201
    };
202 202

	
203 203
    template <typename Graph>
204 204
    struct CountArcsSelector<
205 205
      Graph,
206 206
      typename enable_if<typename Graph::ArcNumTag, void>::type>
207 207
    {
208 208
      static int count(const Graph &g) {
209 209
        return g.arcNum();
210 210
      }
211 211
    };
212 212
  }
213 213

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

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

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

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

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

	
37 34
namespace lemon {
38 35

	
39 36
  /// \addtogroup min_cost_flow
40 37
  /// @{
41 38

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

	
73
    /// The flow type of the algorithm
74
    typedef F Flow;
75
    /// The cost type of the algorithm
75
    /// The type of the flow amounts, capacity bounds and supply values
76
    typedef V Value;
77
    /// The type of the arc costs
76 78
    typedef C Cost;
77
#ifdef DOXYGEN
78
    /// The type of the flow map
79
    typedef GR::ArcMap<Flow> FlowMap;
80
    /// The type of the potential map
81
    typedef GR::NodeMap<Cost> PotentialMap;
82
#else
83
    /// The type of the flow map
84
    typedef typename GR::template ArcMap<Flow> FlowMap;
85
    /// The type of the potential map
86
    typedef typename GR::template NodeMap<Cost> PotentialMap;
87
#endif
88 79

	
89 80
  public:
90 81

	
91
    /// \brief Enum type for selecting the pivot rule.
82
    /// \brief Problem type constants for the \c run() function.
92 83
    ///
93
    /// Enum type for selecting the pivot rule for the \ref run()
84
    /// Enum type containing the problem type constants that can be
85
    /// returned by the \ref run() function of the algorithm.
86
    enum ProblemType {
87
      /// The problem has no feasible solution (flow).
88
      INFEASIBLE,
89
      /// The problem has optimal solution (i.e. it is feasible and
90
      /// bounded), and the algorithm has found optimal flow and node
91
      /// potentials (primal and dual solutions).
92
      OPTIMAL,
93
      /// The objective function of the problem is unbounded, i.e.
94
      /// there is a directed cycle having negative total cost and
95
      /// infinite upper bound.
96
      UNBOUNDED
97
    };
98
    
99
    /// \brief Constants for selecting the type of the supply constraints.
100
    ///
101
    /// Enum type containing constants for selecting the supply type,
102
    /// i.e. the direction of the inequalities in the supply/demand
103
    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
104
    ///
105
    /// The default supply type is \c GEQ, since this form is supported
106
    /// by other minimum cost flow algorithms and the \ref Circulation
107
    /// algorithm, as well.
108
    /// The \c LEQ problem type can be selected using the \ref supplyType()
94 109
    /// function.
95 110
    ///
111
    /// Note that the equality form is a special case of both supply types.
112
    enum SupplyType {
113

	
114
      /// This option means that there are <em>"greater or equal"</em>
115
      /// supply/demand constraints in the definition, i.e. the exact
116
      /// formulation of the problem is the following.
117
      /**
118
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
119
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
120
              sup(u) \quad \forall u\in V \f]
121
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
122
      */
123
      /// It means that the total demand must be greater or equal to the 
124
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
125
      /// negative) and all the supplies have to be carried out from 
126
      /// the supply nodes, but there could be demands that are not 
127
      /// satisfied.
128
      GEQ,
129
      /// It is just an alias for the \c GEQ option.
130
      CARRY_SUPPLIES = GEQ,
131

	
132
      /// This option means that there are <em>"less or equal"</em>
133
      /// supply/demand constraints in the definition, i.e. the exact
134
      /// formulation of the problem is the following.
135
      /**
136
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
137
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
138
              sup(u) \quad \forall u\in V \f]
139
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
140
      */
141
      /// It means that the total demand must be less or equal to the 
142
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
143
      /// positive) and all the demands have to be satisfied, but there
144
      /// could be supplies that are not carried out from the supply
145
      /// nodes.
146
      LEQ,
147
      /// It is just an alias for the \c LEQ option.
148
      SATISFY_DEMANDS = LEQ
149
    };
150
    
151
    /// \brief Constants for selecting the pivot rule.
152
    ///
153
    /// Enum type containing constants for selecting the pivot rule for
154
    /// the \ref run() function.
155
    ///
96 156
    /// \ref NetworkSimplex provides five different pivot rule
97 157
    /// implementations that significantly affect the running time
98 158
    /// of the algorithm.
99 159
    /// By default \ref BLOCK_SEARCH "Block Search" is used, which
100 160
    /// proved to be the most efficient and the most robust on various
101 161
    /// test inputs according to our benchmark tests.
102 162
    /// However another pivot rule can be selected using the \ref run()
103 163
    /// function with the proper parameter.
104 164
    enum PivotRule {
105 165

	
106 166
      /// The First Eligible pivot rule.
107 167
      /// The next eligible arc is selected in a wraparound fashion
108 168
      /// in every iteration.
109 169
      FIRST_ELIGIBLE,
110 170

	
111 171
      /// The Best Eligible pivot rule.
112 172
      /// The best eligible arc is selected in every iteration.
113 173
      BEST_ELIGIBLE,
114 174

	
115 175
      /// The Block Search pivot rule.
116 176
      /// A specified number of arcs are examined in every iteration
117 177
      /// in a wraparound fashion and the best eligible arc is selected
118 178
      /// from this block.
119 179
      BLOCK_SEARCH,
120 180

	
121 181
      /// The Candidate List pivot rule.
122 182
      /// In a major iteration a candidate list is built from eligible arcs
123 183
      /// in a wraparound fashion and in the following minor iterations
124 184
      /// the best eligible arc is selected from this list.
125 185
      CANDIDATE_LIST,
126 186

	
127 187
      /// The Altering Candidate List pivot rule.
128 188
      /// It is a modified version of the Candidate List method.
129 189
      /// It keeps only the several best eligible arcs from the former
130 190
      /// candidate list and extends this list in every iteration.
131 191
      ALTERING_LIST
132 192
    };
133 193
    
134
    /// \brief Enum type for selecting the problem type.
135
    ///
136
    /// Enum type for selecting the problem type, i.e. the direction of
137
    /// the inequalities in the supply/demand constraints of the
138
    /// \ref min_cost_flow "minimum cost flow problem".
139
    ///
140
    /// The default problem type is \c GEQ, since this form is supported
141
    /// by other minimum cost flow algorithms and the \ref Circulation
142
    /// algorithm as well.
143
    /// The \c LEQ problem type can be selected using the \ref problemType()
144
    /// function.
145
    ///
146
    /// Note that the equality form is a special case of both problem type.
147
    enum ProblemType {
148

	
149
      /// This option means that there are "<em>greater or equal</em>"
150
      /// constraints in the defintion, i.e. the exact formulation of the
151
      /// problem is the following.
152
      /**
153
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
154
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
155
              sup(u) \quad \forall u\in V \f]
156
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
157
      */
158
      /// It means that the total demand must be greater or equal to the 
159
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
160
      /// negative) and all the supplies have to be carried out from 
161
      /// the supply nodes, but there could be demands that are not 
162
      /// satisfied.
163
      GEQ,
164
      /// It is just an alias for the \c GEQ option.
165
      CARRY_SUPPLIES = GEQ,
166

	
167
      /// This option means that there are "<em>less or equal</em>"
168
      /// constraints in the defintion, i.e. the exact formulation of the
169
      /// problem is the following.
170
      /**
171
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
172
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
173
              sup(u) \quad \forall u\in V \f]
174
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
175
      */
176
      /// It means that the total demand must be less or equal to the 
177
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
178
      /// positive) and all the demands have to be satisfied, but there
179
      /// could be supplies that are not carried out from the supply
180
      /// nodes.
181
      LEQ,
182
      /// It is just an alias for the \c LEQ option.
183
      SATISFY_DEMANDS = LEQ
184
    };
185

	
186 194
  private:
187 195

	
188 196
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
189 197

	
190
    typedef typename GR::template ArcMap<Flow> FlowArcMap;
191
    typedef typename GR::template ArcMap<Cost> CostArcMap;
192
    typedef typename GR::template NodeMap<Flow> FlowNodeMap;
193

	
194 198
    typedef std::vector<Arc> ArcVector;
195 199
    typedef std::vector<Node> NodeVector;
196 200
    typedef std::vector<int> IntVector;
197 201
    typedef std::vector<bool> BoolVector;
198
    typedef std::vector<Flow> FlowVector;
202
    typedef std::vector<Value> ValueVector;
199 203
    typedef std::vector<Cost> CostVector;
200 204

	
201 205
    // State constants for arcs
202 206
    enum ArcStateEnum {
203 207
      STATE_UPPER = -1,
204 208
      STATE_TREE  =  0,
205 209
      STATE_LOWER =  1
206 210
    };
207 211

	
208 212
  private:
209 213

	
210 214
    // Data related to the underlying digraph
211 215
    const GR &_graph;
212 216
    int _node_num;
213 217
    int _arc_num;
214 218

	
215 219
    // Parameters of the problem
216
    FlowArcMap *_plower;
217
    FlowArcMap *_pupper;
218
    CostArcMap *_pcost;
219
    FlowNodeMap *_psupply;
220
    bool _pstsup;
221
    Node _psource, _ptarget;
222
    Flow _pstflow;
223
    ProblemType _ptype;
224

	
225
    // Result maps
226
    FlowMap *_flow_map;
227
    PotentialMap *_potential_map;
228
    bool _local_flow;
229
    bool _local_potential;
220
    bool _have_lower;
221
    SupplyType _stype;
222
    Value _sum_supply;
230 223

	
231 224
    // Data structures for storing the digraph
232 225
    IntNodeMap _node_id;
233
    ArcVector _arc_ref;
226
    IntArcMap _arc_id;
234 227
    IntVector _source;
235 228
    IntVector _target;
236 229

	
237 230
    // Node and arc data
238
    FlowVector _cap;
231
    ValueVector _lower;
232
    ValueVector _upper;
233
    ValueVector _cap;
239 234
    CostVector _cost;
240
    FlowVector _supply;
241
    FlowVector _flow;
235
    ValueVector _supply;
236
    ValueVector _flow;
242 237
    CostVector _pi;
243 238

	
244 239
    // Data for storing the spanning tree structure
245 240
    IntVector _parent;
246 241
    IntVector _pred;
247 242
    IntVector _thread;
248 243
    IntVector _rev_thread;
249 244
    IntVector _succ_num;
250 245
    IntVector _last_succ;
251 246
    IntVector _dirty_revs;
252 247
    BoolVector _forward;
253 248
    IntVector _state;
254 249
    int _root;
255 250

	
256 251
    // Temporary data used in the current pivot iteration
257 252
    int in_arc, join, u_in, v_in, u_out, v_out;
258 253
    int first, second, right, last;
259 254
    int stem, par_stem, new_stem;
260
    Flow delta;
255
    Value delta;
256

	
257
  public:
258
  
259
    /// \brief Constant for infinite upper bounds (capacities).
260
    ///
261
    /// Constant for infinite upper bounds (capacities).
262
    /// It is \c std::numeric_limits<Value>::infinity() if available,
263
    /// \c std::numeric_limits<Value>::max() otherwise.
264
    const Value INF;
261 265

	
262 266
  private:
263 267

	
264 268
    // Implementation of the First Eligible pivot rule
265 269
    class FirstEligiblePivotRule
266 270
    {
267 271
    private:
268 272

	
269 273
      // References to the NetworkSimplex class
270 274
      const IntVector  &_source;
271 275
      const IntVector  &_target;
272 276
      const CostVector &_cost;
273 277
      const IntVector  &_state;
274 278
      const CostVector &_pi;
275 279
      int &_in_arc;
276 280
      int _arc_num;
277 281

	
278 282
      // Pivot rule data
279 283
      int _next_arc;
280 284

	
281 285
    public:
282 286

	
283 287
      // Constructor
284 288
      FirstEligiblePivotRule(NetworkSimplex &ns) :
285 289
        _source(ns._source), _target(ns._target),
286 290
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
287 291
        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
288 292
      {}
289 293

	
290 294
      // Find next entering arc
291 295
      bool findEnteringArc() {
292 296
        Cost c;
293 297
        for (int e = _next_arc; e < _arc_num; ++e) {
294 298
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
295 299
          if (c < 0) {
296 300
            _in_arc = e;
297 301
            _next_arc = e + 1;
298 302
            return true;
299 303
          }
300 304
        }
301 305
        for (int e = 0; e < _next_arc; ++e) {
302 306
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
303 307
          if (c < 0) {
304 308
            _in_arc = e;
305 309
            _next_arc = e + 1;
306 310
            return true;
307 311
          }
308 312
        }
309 313
        return false;
310 314
      }
311 315

	
312 316
    }; //class FirstEligiblePivotRule
313 317

	
314 318

	
315 319
    // Implementation of the Best Eligible pivot rule
316 320
    class BestEligiblePivotRule
317 321
    {
318 322
    private:
319 323

	
320 324
      // References to the NetworkSimplex class
321 325
      const IntVector  &_source;
322 326
      const IntVector  &_target;
323 327
      const CostVector &_cost;
324 328
      const IntVector  &_state;
325 329
      const CostVector &_pi;
326 330
      int &_in_arc;
327 331
      int _arc_num;
328 332

	
329 333
    public:
330 334

	
331 335
      // Constructor
332 336
      BestEligiblePivotRule(NetworkSimplex &ns) :
333 337
        _source(ns._source), _target(ns._target),
334 338
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
335 339
        _in_arc(ns.in_arc), _arc_num(ns._arc_num)
336 340
      {}
337 341

	
338 342
      // Find next entering arc
339 343
      bool findEnteringArc() {
340 344
        Cost c, min = 0;
341 345
        for (int e = 0; e < _arc_num; ++e) {
342 346
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
343 347
          if (c < min) {
344 348
            min = c;
345 349
            _in_arc = e;
346 350
          }
347 351
        }
348 352
        return min < 0;
349 353
      }
350 354

	
351 355
    }; //class BestEligiblePivotRule
352 356

	
353 357

	
354 358
    // Implementation of the Block Search pivot rule
355 359
    class BlockSearchPivotRule
356 360
    {
357 361
    private:
358 362

	
359 363
      // References to the NetworkSimplex class
360 364
      const IntVector  &_source;
361 365
      const IntVector  &_target;
362 366
      const CostVector &_cost;
363 367
      const IntVector  &_state;
364 368
      const CostVector &_pi;
365 369
      int &_in_arc;
366 370
      int _arc_num;
367 371

	
368 372
      // Pivot rule data
369 373
      int _block_size;
370 374
      int _next_arc;
371 375

	
372 376
    public:
373 377

	
374 378
      // Constructor
375 379
      BlockSearchPivotRule(NetworkSimplex &ns) :
376 380
        _source(ns._source), _target(ns._target),
377 381
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
378 382
        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
379 383
      {
380 384
        // The main parameters of the pivot rule
381 385
        const double BLOCK_SIZE_FACTOR = 2.0;
382 386
        const int MIN_BLOCK_SIZE = 10;
383 387

	
384 388
        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
385 389
                                    std::sqrt(double(_arc_num))),
386 390
                                MIN_BLOCK_SIZE );
387 391
      }
388 392

	
389 393
      // Find next entering arc
390 394
      bool findEnteringArc() {
391 395
        Cost c, min = 0;
392 396
        int cnt = _block_size;
393 397
        int e, min_arc = _next_arc;
394 398
        for (e = _next_arc; e < _arc_num; ++e) {
395 399
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
396 400
          if (c < min) {
397 401
            min = c;
398 402
            min_arc = e;
399 403
          }
400 404
          if (--cnt == 0) {
401 405
            if (min < 0) break;
402 406
            cnt = _block_size;
403 407
          }
404 408
        }
405 409
        if (min == 0 || cnt > 0) {
406 410
          for (e = 0; e < _next_arc; ++e) {
407 411
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
408 412
            if (c < min) {
409 413
              min = c;
410 414
              min_arc = e;
411 415
            }
412 416
            if (--cnt == 0) {
413 417
              if (min < 0) break;
414 418
              cnt = _block_size;
415 419
            }
416 420
          }
417 421
        }
418 422
        if (min >= 0) return false;
419 423
        _in_arc = min_arc;
420 424
        _next_arc = e;
421 425
        return true;
422 426
      }
423 427

	
424 428
    }; //class BlockSearchPivotRule
425 429

	
426 430

	
427 431
    // Implementation of the Candidate List pivot rule
428 432
    class CandidateListPivotRule
429 433
    {
430 434
    private:
431 435

	
432 436
      // References to the NetworkSimplex class
433 437
      const IntVector  &_source;
434 438
      const IntVector  &_target;
435 439
      const CostVector &_cost;
436 440
      const IntVector  &_state;
437 441
      const CostVector &_pi;
438 442
      int &_in_arc;
439 443
      int _arc_num;
440 444

	
441 445
      // Pivot rule data
442 446
      IntVector _candidates;
443 447
      int _list_length, _minor_limit;
444 448
      int _curr_length, _minor_count;
445 449
      int _next_arc;
446 450

	
447 451
    public:
448 452

	
449 453
      /// Constructor
450 454
      CandidateListPivotRule(NetworkSimplex &ns) :
451 455
        _source(ns._source), _target(ns._target),
452 456
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
... ...
@@ -470,1117 +474,941 @@
470 474
      /// Find next entering arc
471 475
      bool findEnteringArc() {
472 476
        Cost min, c;
473 477
        int e, min_arc = _next_arc;
474 478
        if (_curr_length > 0 && _minor_count < _minor_limit) {
475 479
          // Minor iteration: select the best eligible arc from the
476 480
          // current candidate list
477 481
          ++_minor_count;
478 482
          min = 0;
479 483
          for (int i = 0; i < _curr_length; ++i) {
480 484
            e = _candidates[i];
481 485
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
482 486
            if (c < min) {
483 487
              min = c;
484 488
              min_arc = e;
485 489
            }
486 490
            if (c >= 0) {
487 491
              _candidates[i--] = _candidates[--_curr_length];
488 492
            }
489 493
          }
490 494
          if (min < 0) {
491 495
            _in_arc = min_arc;
492 496
            return true;
493 497
          }
494 498
        }
495 499

	
496 500
        // Major iteration: build a new candidate list
497 501
        min = 0;
498 502
        _curr_length = 0;
499 503
        for (e = _next_arc; e < _arc_num; ++e) {
500 504
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
501 505
          if (c < 0) {
502 506
            _candidates[_curr_length++] = e;
503 507
            if (c < min) {
504 508
              min = c;
505 509
              min_arc = e;
506 510
            }
507 511
            if (_curr_length == _list_length) break;
508 512
          }
509 513
        }
510 514
        if (_curr_length < _list_length) {
511 515
          for (e = 0; e < _next_arc; ++e) {
512 516
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
513 517
            if (c < 0) {
514 518
              _candidates[_curr_length++] = e;
515 519
              if (c < min) {
516 520
                min = c;
517 521
                min_arc = e;
518 522
              }
519 523
              if (_curr_length == _list_length) break;
520 524
            }
521 525
          }
522 526
        }
523 527
        if (_curr_length == 0) return false;
524 528
        _minor_count = 1;
525 529
        _in_arc = min_arc;
526 530
        _next_arc = e;
527 531
        return true;
528 532
      }
529 533

	
530 534
    }; //class CandidateListPivotRule
531 535

	
532 536

	
533 537
    // Implementation of the Altering Candidate List pivot rule
534 538
    class AlteringListPivotRule
535 539
    {
536 540
    private:
537 541

	
538 542
      // References to the NetworkSimplex class
539 543
      const IntVector  &_source;
540 544
      const IntVector  &_target;
541 545
      const CostVector &_cost;
542 546
      const IntVector  &_state;
543 547
      const CostVector &_pi;
544 548
      int &_in_arc;
545 549
      int _arc_num;
546 550

	
547 551
      // Pivot rule data
548 552
      int _block_size, _head_length, _curr_length;
549 553
      int _next_arc;
550 554
      IntVector _candidates;
551 555
      CostVector _cand_cost;
552 556

	
553 557
      // Functor class to compare arcs during sort of the candidate list
554 558
      class SortFunc
555 559
      {
556 560
      private:
557 561
        const CostVector &_map;
558 562
      public:
559 563
        SortFunc(const CostVector &map) : _map(map) {}
560 564
        bool operator()(int left, int right) {
561 565
          return _map[left] > _map[right];
562 566
        }
563 567
      };
564 568

	
565 569
      SortFunc _sort_func;
566 570

	
567 571
    public:
568 572

	
569 573
      // Constructor
570 574
      AlteringListPivotRule(NetworkSimplex &ns) :
571 575
        _source(ns._source), _target(ns._target),
572 576
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
573 577
        _in_arc(ns.in_arc), _arc_num(ns._arc_num),
574 578
        _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
575 579
      {
576 580
        // The main parameters of the pivot rule
577 581
        const double BLOCK_SIZE_FACTOR = 1.5;
578 582
        const int MIN_BLOCK_SIZE = 10;
579 583
        const double HEAD_LENGTH_FACTOR = 0.1;
580 584
        const int MIN_HEAD_LENGTH = 3;
581 585

	
582 586
        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
583 587
                                    std::sqrt(double(_arc_num))),
584 588
                                MIN_BLOCK_SIZE );
585 589
        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
586 590
                                 MIN_HEAD_LENGTH );
587 591
        _candidates.resize(_head_length + _block_size);
588 592
        _curr_length = 0;
589 593
      }
590 594

	
591 595
      // Find next entering arc
592 596
      bool findEnteringArc() {
593 597
        // Check the current candidate list
594 598
        int e;
595 599
        for (int i = 0; i < _curr_length; ++i) {
596 600
          e = _candidates[i];
597 601
          _cand_cost[e] = _state[e] *
598 602
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
599 603
          if (_cand_cost[e] >= 0) {
600 604
            _candidates[i--] = _candidates[--_curr_length];
601 605
          }
602 606
        }
603 607

	
604 608
        // Extend the list
605 609
        int cnt = _block_size;
606 610
        int last_arc = 0;
607 611
        int limit = _head_length;
608 612

	
609 613
        for (int e = _next_arc; e < _arc_num; ++e) {
610 614
          _cand_cost[e] = _state[e] *
611 615
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
612 616
          if (_cand_cost[e] < 0) {
613 617
            _candidates[_curr_length++] = e;
614 618
            last_arc = e;
615 619
          }
616 620
          if (--cnt == 0) {
617 621
            if (_curr_length > limit) break;
618 622
            limit = 0;
619 623
            cnt = _block_size;
620 624
          }
621 625
        }
622 626
        if (_curr_length <= limit) {
623 627
          for (int e = 0; e < _next_arc; ++e) {
624 628
            _cand_cost[e] = _state[e] *
625 629
              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
626 630
            if (_cand_cost[e] < 0) {
627 631
              _candidates[_curr_length++] = e;
628 632
              last_arc = e;
629 633
            }
630 634
            if (--cnt == 0) {
631 635
              if (_curr_length > limit) break;
632 636
              limit = 0;
633 637
              cnt = _block_size;
634 638
            }
635 639
          }
636 640
        }
637 641
        if (_curr_length == 0) return false;
638 642
        _next_arc = last_arc + 1;
639 643

	
640 644
        // Make heap of the candidate list (approximating a partial sort)
641 645
        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
642 646
                   _sort_func );
643 647

	
644 648
        // Pop the first element of the heap
645 649
        _in_arc = _candidates[0];
646 650
        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
647 651
                  _sort_func );
648 652
        _curr_length = std::min(_head_length, _curr_length - 1);
649 653
        return true;
650 654
      }
651 655

	
652 656
    }; //class AlteringListPivotRule
653 657

	
654 658
  public:
655 659

	
656 660
    /// \brief Constructor.
657 661
    ///
658 662
    /// The constructor of the class.
659 663
    ///
660 664
    /// \param graph The digraph the algorithm runs on.
661 665
    NetworkSimplex(const GR& graph) :
662
      _graph(graph),
663
      _plower(NULL), _pupper(NULL), _pcost(NULL),
664
      _psupply(NULL), _pstsup(false), _ptype(GEQ),
665
      _flow_map(NULL), _potential_map(NULL),
666
      _local_flow(false), _local_potential(false),
667
      _node_id(graph)
666
      _graph(graph), _node_id(graph), _arc_id(graph),
667
      INF(std::numeric_limits<Value>::has_infinity ?
668
          std::numeric_limits<Value>::infinity() :
669
          std::numeric_limits<Value>::max())
668 670
    {
669
      LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
670
                   std::numeric_limits<Flow>::is_signed,
671
        "The flow type of NetworkSimplex must be signed integer");
672
      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
673
                   std::numeric_limits<Cost>::is_signed,
674
        "The cost type of NetworkSimplex must be signed integer");
675
    }
671
      // Check the value types
672
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
673
        "The flow type of NetworkSimplex must be signed");
674
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
675
        "The cost type of NetworkSimplex must be signed");
676
        
677
      // Resize vectors
678
      _node_num = countNodes(_graph);
679
      _arc_num = countArcs(_graph);
680
      int all_node_num = _node_num + 1;
681
      int all_arc_num = _arc_num + _node_num;
676 682

	
677
    /// Destructor.
678
    ~NetworkSimplex() {
679
      if (_local_flow) delete _flow_map;
680
      if (_local_potential) delete _potential_map;
683
      _source.resize(all_arc_num);
684
      _target.resize(all_arc_num);
685

	
686
      _lower.resize(all_arc_num);
687
      _upper.resize(all_arc_num);
688
      _cap.resize(all_arc_num);
689
      _cost.resize(all_arc_num);
690
      _supply.resize(all_node_num);
691
      _flow.resize(all_arc_num);
692
      _pi.resize(all_node_num);
693

	
694
      _parent.resize(all_node_num);
695
      _pred.resize(all_node_num);
696
      _forward.resize(all_node_num);
697
      _thread.resize(all_node_num);
698
      _rev_thread.resize(all_node_num);
699
      _succ_num.resize(all_node_num);
700
      _last_succ.resize(all_node_num);
701
      _state.resize(all_arc_num);
702

	
703
      // Copy the graph (store the arcs in a mixed order)
704
      int i = 0;
705
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
706
        _node_id[n] = i;
707
      }
708
      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
709
      i = 0;
710
      for (ArcIt a(_graph); a != INVALID; ++a) {
711
        _arc_id[a] = i;
712
        _source[i] = _node_id[_graph.source(a)];
713
        _target[i] = _node_id[_graph.target(a)];
714
        if ((i += k) >= _arc_num) i = (i % k) + 1;
715
      }
716
      
717
      // Initialize maps
718
      for (int i = 0; i != _node_num; ++i) {
719
        _supply[i] = 0;
720
      }
721
      for (int i = 0; i != _arc_num; ++i) {
722
        _lower[i] = 0;
723
        _upper[i] = INF;
724
        _cost[i] = 1;
725
      }
726
      _have_lower = false;
727
      _stype = GEQ;
681 728
    }
682 729

	
683 730
    /// \name Parameters
684 731
    /// The parameters of the algorithm can be specified using these
685 732
    /// functions.
686 733

	
687 734
    /// @{
688 735

	
689 736
    /// \brief Set the lower bounds on the arcs.
690 737
    ///
691 738
    /// This function sets the lower bounds on the arcs.
692
    /// If neither this function nor \ref boundMaps() is used before
693
    /// calling \ref run(), the lower bounds will be set to zero
694
    /// on all arcs.
739
    /// If it is not used before calling \ref run(), the lower bounds
740
    /// will be set to zero on all arcs.
695 741
    ///
696 742
    /// \param map An arc map storing the lower bounds.
697
    /// Its \c Value type must be convertible to the \c Flow type
743
    /// Its \c Value type must be convertible to the \c Value type
698 744
    /// of the algorithm.
699 745
    ///
700 746
    /// \return <tt>(*this)</tt>
701
    template <typename LOWER>
702
    NetworkSimplex& lowerMap(const LOWER& map) {
703
      delete _plower;
704
      _plower = new FlowArcMap(_graph);
747
    template <typename LowerMap>
748
    NetworkSimplex& lowerMap(const LowerMap& map) {
749
      _have_lower = true;
705 750
      for (ArcIt a(_graph); a != INVALID; ++a) {
706
        (*_plower)[a] = map[a];
751
        _lower[_arc_id[a]] = map[a];
707 752
      }
708 753
      return *this;
709 754
    }
710 755

	
711 756
    /// \brief Set the upper bounds (capacities) on the arcs.
712 757
    ///
713 758
    /// This function sets the upper bounds (capacities) on the arcs.
714
    /// If none of the functions \ref upperMap(), \ref capacityMap()
715
    /// and \ref boundMaps() is used before calling \ref run(),
716
    /// the upper bounds (capacities) will be set to
717
    /// \c std::numeric_limits<Flow>::max() on all arcs.
759
    /// If it is not used before calling \ref run(), the upper bounds
760
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
761
    /// unbounded from above on each arc).
718 762
    ///
719 763
    /// \param map An arc map storing the upper bounds.
720
    /// Its \c Value type must be convertible to the \c Flow type
764
    /// Its \c Value type must be convertible to the \c Value type
721 765
    /// of the algorithm.
722 766
    ///
723 767
    /// \return <tt>(*this)</tt>
724
    template<typename UPPER>
725
    NetworkSimplex& upperMap(const UPPER& map) {
726
      delete _pupper;
727
      _pupper = new FlowArcMap(_graph);
768
    template<typename UpperMap>
769
    NetworkSimplex& upperMap(const UpperMap& map) {
728 770
      for (ArcIt a(_graph); a != INVALID; ++a) {
729
        (*_pupper)[a] = map[a];
771
        _upper[_arc_id[a]] = map[a];
730 772
      }
731 773
      return *this;
732 774
    }
733 775

	
734
    /// \brief Set the upper bounds (capacities) on the arcs.
735
    ///
736
    /// This function sets the upper bounds (capacities) on the arcs.
737
    /// It is just an alias for \ref upperMap().
738
    ///
739
    /// \return <tt>(*this)</tt>
740
    template<typename CAP>
741
    NetworkSimplex& capacityMap(const CAP& map) {
742
      return upperMap(map);
743
    }
744

	
745
    /// \brief Set the lower and upper bounds on the arcs.
746
    ///
747
    /// This function sets the lower and upper bounds on the arcs.
748
    /// If neither this function nor \ref lowerMap() is used before
749
    /// calling \ref run(), the lower bounds will be set to zero
750
    /// on all arcs.
751
    /// If none of the functions \ref upperMap(), \ref capacityMap()
752
    /// and \ref boundMaps() is used before calling \ref run(),
753
    /// the upper bounds (capacities) will be set to
754
    /// \c std::numeric_limits<Flow>::max() on all arcs.
755
    ///
756
    /// \param lower An arc map storing the lower bounds.
757
    /// \param upper An arc map storing the upper bounds.
758
    ///
759
    /// The \c Value type of the maps must be convertible to the
760
    /// \c Flow type of the algorithm.
761
    ///
762
    /// \note This function is just a shortcut of calling \ref lowerMap()
763
    /// and \ref upperMap() separately.
764
    ///
765
    /// \return <tt>(*this)</tt>
766
    template <typename LOWER, typename UPPER>
767
    NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
768
      return lowerMap(lower).upperMap(upper);
769
    }
770

	
771 776
    /// \brief Set the costs of the arcs.
772 777
    ///
773 778
    /// This function sets the costs of the arcs.
774 779
    /// If it is not used before calling \ref run(), the costs
775 780
    /// will be set to \c 1 on all arcs.
776 781
    ///
777 782
    /// \param map An arc map storing the costs.
778 783
    /// Its \c Value type must be convertible to the \c Cost type
779 784
    /// of the algorithm.
780 785
    ///
781 786
    /// \return <tt>(*this)</tt>
782
    template<typename COST>
783
    NetworkSimplex& costMap(const COST& map) {
784
      delete _pcost;
785
      _pcost = new CostArcMap(_graph);
787
    template<typename CostMap>
788
    NetworkSimplex& costMap(const CostMap& map) {
786 789
      for (ArcIt a(_graph); a != INVALID; ++a) {
787
        (*_pcost)[a] = map[a];
790
        _cost[_arc_id[a]] = map[a];
788 791
      }
789 792
      return *this;
790 793
    }
791 794

	
792 795
    /// \brief Set the supply values of the nodes.
793 796
    ///
794 797
    /// This function sets the supply values of the nodes.
795 798
    /// If neither this function nor \ref stSupply() is used before
796 799
    /// calling \ref run(), the supply of each node will be set to zero.
797 800
    /// (It makes sense only if non-zero lower bounds are given.)
798 801
    ///
799 802
    /// \param map A node map storing the supply values.
800
    /// Its \c Value type must be convertible to the \c Flow type
803
    /// Its \c Value type must be convertible to the \c Value type
801 804
    /// of the algorithm.
802 805
    ///
803 806
    /// \return <tt>(*this)</tt>
804
    template<typename SUP>
805
    NetworkSimplex& supplyMap(const SUP& map) {
806
      delete _psupply;
807
      _pstsup = false;
808
      _psupply = new FlowNodeMap(_graph);
807
    template<typename SupplyMap>
808
    NetworkSimplex& supplyMap(const SupplyMap& map) {
809 809
      for (NodeIt n(_graph); n != INVALID; ++n) {
810
        (*_psupply)[n] = map[n];
810
        _supply[_node_id[n]] = map[n];
811 811
      }
812 812
      return *this;
813 813
    }
814 814

	
815 815
    /// \brief Set single source and target nodes and a supply value.
816 816
    ///
817 817
    /// This function sets a single source node and a single target node
818 818
    /// and the required flow value.
819 819
    /// If neither this function nor \ref supplyMap() is used before
820 820
    /// calling \ref run(), the supply of each node will be set to zero.
821 821
    /// (It makes sense only if non-zero lower bounds are given.)
822 822
    ///
823
    /// Using this function has the same effect as using \ref supplyMap()
824
    /// with such a map in which \c k is assigned to \c s, \c -k is
825
    /// assigned to \c t and all other nodes have zero supply value.
826
    ///
823 827
    /// \param s The source node.
824 828
    /// \param t The target node.
825 829
    /// \param k The required amount of flow from node \c s to node \c t
826 830
    /// (i.e. the supply of \c s and the demand of \c t).
827 831
    ///
828 832
    /// \return <tt>(*this)</tt>
829
    NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
830
      delete _psupply;
831
      _psupply = NULL;
832
      _pstsup = true;
833
      _psource = s;
834
      _ptarget = t;
835
      _pstflow = k;
833
    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
834
      for (int i = 0; i != _node_num; ++i) {
835
        _supply[i] = 0;
836
      }
837
      _supply[_node_id[s]] =  k;
838
      _supply[_node_id[t]] = -k;
836 839
      return *this;
837 840
    }
838 841
    
839
    /// \brief Set the problem type.
842
    /// \brief Set the type of the supply constraints.
840 843
    ///
841
    /// This function sets the problem type for the algorithm.
842
    /// If it is not used before calling \ref run(), the \ref GEQ problem
844
    /// This function sets the type of the supply/demand constraints.
845
    /// If it is not used before calling \ref run(), the \ref GEQ supply
843 846
    /// type will be used.
844 847
    ///
845
    /// For more information see \ref ProblemType.
848
    /// For more information see \ref SupplyType.
846 849
    ///
847 850
    /// \return <tt>(*this)</tt>
848
    NetworkSimplex& problemType(ProblemType problem_type) {
849
      _ptype = problem_type;
851
    NetworkSimplex& supplyType(SupplyType supply_type) {
852
      _stype = supply_type;
850 853
      return *this;
851 854
    }
852 855

	
853
    /// \brief Set the flow map.
854
    ///
855
    /// This function sets the flow map.
856
    /// If it is not used before calling \ref run(), an instance will
857
    /// be allocated automatically. The destructor deallocates this
858
    /// automatically allocated map, of course.
859
    ///
860
    /// \return <tt>(*this)</tt>
861
    NetworkSimplex& flowMap(FlowMap& map) {
862
      if (_local_flow) {
863
        delete _flow_map;
864
        _local_flow = false;
865
      }
866
      _flow_map = &map;
867
      return *this;
868
    }
869

	
870
    /// \brief Set the potential map.
871
    ///
872
    /// This function sets the potential map, which is used for storing
873
    /// the dual solution.
874
    /// If it is not used before calling \ref run(), an instance will
875
    /// be allocated automatically. The destructor deallocates this
876
    /// automatically allocated map, of course.
877
    ///
878
    /// \return <tt>(*this)</tt>
879
    NetworkSimplex& potentialMap(PotentialMap& map) {
880
      if (_local_potential) {
881
        delete _potential_map;
882
        _local_potential = false;
883
      }
884
      _potential_map = &map;
885
      return *this;
886
    }
887
    
888 856
    /// @}
889 857

	
890 858
    /// \name Execution Control
891 859
    /// The algorithm can be executed using \ref run().
892 860

	
893 861
    /// @{
894 862

	
895 863
    /// \brief Run the algorithm.
896 864
    ///
897 865
    /// This function runs the algorithm.
898 866
    /// The paramters can be specified using functions \ref lowerMap(),
899
    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
900
    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
901
    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
867
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 
868
    /// \ref supplyType().
902 869
    /// For example,
903 870
    /// \code
904 871
    ///   NetworkSimplex<ListDigraph> ns(graph);
905
    ///   ns.boundMaps(lower, upper).costMap(cost)
872
    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
906 873
    ///     .supplyMap(sup).run();
907 874
    /// \endcode
908 875
    ///
909 876
    /// This function can be called more than once. All the parameters
910 877
    /// that have been given are kept for the next call, unless
911 878
    /// \ref reset() is called, thus only the modified parameters
912 879
    /// have to be set again. See \ref reset() for examples.
880
    /// However the underlying digraph must not be modified after this
881
    /// class have been constructed, since it copies and extends the graph.
913 882
    ///
914 883
    /// \param pivot_rule The pivot rule that will be used during the
915 884
    /// algorithm. For more information see \ref PivotRule.
916 885
    ///
917
    /// \return \c true if a feasible flow can be found.
918
    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
919
      return init() && start(pivot_rule);
886
    /// \return \c INFEASIBLE if no feasible flow exists,
887
    /// \n \c OPTIMAL if the problem has optimal solution
888
    /// (i.e. it is feasible and bounded), and the algorithm has found
889
    /// optimal flow and node potentials (primal and dual solutions),
890
    /// \n \c UNBOUNDED if the objective function of the problem is
891
    /// unbounded, i.e. there is a directed cycle having negative total
892
    /// cost and infinite upper bound.
893
    ///
894
    /// \see ProblemType, PivotRule
895
    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
896
      if (!init()) return INFEASIBLE;
897
      return start(pivot_rule);
920 898
    }
921 899

	
922 900
    /// \brief Reset all the parameters that have been given before.
923 901
    ///
924 902
    /// This function resets all the paramaters that have been given
925 903
    /// before using functions \ref lowerMap(), \ref upperMap(),
926
    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
927
    /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
928
    /// \ref flowMap() and \ref potentialMap().
904
    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
929 905
    ///
930 906
    /// It is useful for multiple run() calls. If this function is not
931 907
    /// used, all the parameters given before are kept for the next
932 908
    /// \ref run() call.
909
    /// However the underlying digraph must not be modified after this
910
    /// class have been constructed, since it copies and extends the graph.
933 911
    ///
934 912
    /// For example,
935 913
    /// \code
936 914
    ///   NetworkSimplex<ListDigraph> ns(graph);
937 915
    ///
938 916
    ///   // First run
939
    ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
917
    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
940 918
    ///     .supplyMap(sup).run();
941 919
    ///
942 920
    ///   // Run again with modified cost map (reset() is not called,
943 921
    ///   // so only the cost map have to be set again)
944 922
    ///   cost[e] += 100;
945 923
    ///   ns.costMap(cost).run();
946 924
    ///
947 925
    ///   // Run again from scratch using reset()
948 926
    ///   // (the lower bounds will be set to zero on all arcs)
949 927
    ///   ns.reset();
950
    ///   ns.capacityMap(cap).costMap(cost)
928
    ///   ns.upperMap(capacity).costMap(cost)
951 929
    ///     .supplyMap(sup).run();
952 930
    /// \endcode
953 931
    ///
954 932
    /// \return <tt>(*this)</tt>
955 933
    NetworkSimplex& reset() {
956
      delete _plower;
957
      delete _pupper;
958
      delete _pcost;
959
      delete _psupply;
960
      _plower = NULL;
961
      _pupper = NULL;
962
      _pcost = NULL;
963
      _psupply = NULL;
964
      _pstsup = false;
965
      _ptype = GEQ;
966
      if (_local_flow) delete _flow_map;
967
      if (_local_potential) delete _potential_map;
968
      _flow_map = NULL;
969
      _potential_map = NULL;
970
      _local_flow = false;
971
      _local_potential = false;
972

	
934
      for (int i = 0; i != _node_num; ++i) {
935
        _supply[i] = 0;
936
      }
937
      for (int i = 0; i != _arc_num; ++i) {
938
        _lower[i] = 0;
939
        _upper[i] = INF;
940
        _cost[i] = 1;
941
      }
942
      _have_lower = false;
943
      _stype = GEQ;
973 944
      return *this;
974 945
    }
975 946

	
976 947
    /// @}
977 948

	
978 949
    /// \name Query Functions
979 950
    /// The results of the algorithm can be obtained using these
980 951
    /// functions.\n
981 952
    /// The \ref run() function must be called before using them.
982 953

	
983 954
    /// @{
984 955

	
985 956
    /// \brief Return the total cost of the found flow.
986 957
    ///
987 958
    /// This function returns the total cost of the found flow.
988
    /// The complexity of the function is O(e).
959
    /// Its complexity is O(e).
989 960
    ///
990 961
    /// \note The return type of the function can be specified as a
991 962
    /// template parameter. For example,
992 963
    /// \code
993 964
    ///   ns.totalCost<double>();
994 965
    /// \endcode
995 966
    /// It is useful if the total cost cannot be stored in the \c Cost
996 967
    /// type of the algorithm, which is the default return type of the
997 968
    /// function.
998 969
    ///
999 970
    /// \pre \ref run() must be called before using this function.
1000
    template <typename Num>
1001
    Num totalCost() const {
1002
      Num c = 0;
1003
      if (_pcost) {
1004
        for (ArcIt e(_graph); e != INVALID; ++e)
1005
          c += (*_flow_map)[e] * (*_pcost)[e];
1006
      } else {
1007
        for (ArcIt e(_graph); e != INVALID; ++e)
1008
          c += (*_flow_map)[e];
971
    template <typename Number>
972
    Number totalCost() const {
973
      Number c = 0;
974
      for (ArcIt a(_graph); a != INVALID; ++a) {
975
        int i = _arc_id[a];
976
        c += Number(_flow[i]) * Number(_cost[i]);
1009 977
      }
1010 978
      return c;
1011 979
    }
1012 980

	
1013 981
#ifndef DOXYGEN
1014 982
    Cost totalCost() const {
1015 983
      return totalCost<Cost>();
1016 984
    }
1017 985
#endif
1018 986

	
1019 987
    /// \brief Return the flow on the given arc.
1020 988
    ///
1021 989
    /// This function returns the flow on the given arc.
1022 990
    ///
1023 991
    /// \pre \ref run() must be called before using this function.
1024
    Flow flow(const Arc& a) const {
1025
      return (*_flow_map)[a];
992
    Value flow(const Arc& a) const {
993
      return _flow[_arc_id[a]];
1026 994
    }
1027 995

	
1028
    /// \brief Return a const reference to the flow map.
996
    /// \brief Return the flow map (the primal solution).
1029 997
    ///
1030
    /// This function returns a const reference to an arc map storing
1031
    /// the found flow.
998
    /// This function copies the flow value on each arc into the given
999
    /// map. The \c Value type of the algorithm must be convertible to
1000
    /// the \c Value type of the map.
1032 1001
    ///
1033 1002
    /// \pre \ref run() must be called before using this function.
1034
    const FlowMap& flowMap() const {
1035
      return *_flow_map;
1003
    template <typename FlowMap>
1004
    void flowMap(FlowMap &map) const {
1005
      for (ArcIt a(_graph); a != INVALID; ++a) {
1006
        map.set(a, _flow[_arc_id[a]]);
1007
      }
1036 1008
    }
1037 1009

	
1038 1010
    /// \brief Return the potential (dual value) of the given node.
1039 1011
    ///
1040 1012
    /// This function returns the potential (dual value) of the
1041 1013
    /// given node.
1042 1014
    ///
1043 1015
    /// \pre \ref run() must be called before using this function.
1044 1016
    Cost potential(const Node& n) const {
1045
      return (*_potential_map)[n];
1017
      return _pi[_node_id[n]];
1046 1018
    }
1047 1019

	
1048
    /// \brief Return a const reference to the potential map
1049
    /// (the dual solution).
1020
    /// \brief Return the potential map (the dual solution).
1050 1021
    ///
1051
    /// This function returns a const reference to a node map storing
1052
    /// the found potentials, which form the dual solution of the
1053
    /// \ref min_cost_flow "minimum cost flow" problem.
1022
    /// This function copies the potential (dual value) of each node
1023
    /// into the given map.
1024
    /// The \c Cost type of the algorithm must be convertible to the
1025
    /// \c Value type of the map.
1054 1026
    ///
1055 1027
    /// \pre \ref run() must be called before using this function.
1056
    const PotentialMap& potentialMap() const {
1057
      return *_potential_map;
1028
    template <typename PotentialMap>
1029
    void potentialMap(PotentialMap &map) const {
1030
      for (NodeIt n(_graph); n != INVALID; ++n) {
1031
        map.set(n, _pi[_node_id[n]]);
1032
      }
1058 1033
    }
1059 1034

	
1060 1035
    /// @}
1061 1036

	
1062 1037
  private:
1063 1038

	
1064 1039
    // Initialize internal data structures
1065 1040
    bool init() {
1066
      // Initialize result maps
1067
      if (!_flow_map) {
1068
        _flow_map = new FlowMap(_graph);
1069
        _local_flow = true;
1041
      if (_node_num == 0) return false;
1042

	
1043
      // Check the sum of supply values
1044
      _sum_supply = 0;
1045
      for (int i = 0; i != _node_num; ++i) {
1046
        _sum_supply += _supply[i];
1070 1047
      }
1071
      if (!_potential_map) {
1072
        _potential_map = new PotentialMap(_graph);
1073
        _local_potential = true;
1048
      if ( !((_stype == GEQ && _sum_supply <= 0) ||
1049
             (_stype == LEQ && _sum_supply >= 0)) ) return false;
1050

	
1051
      // Remove non-zero lower bounds
1052
      if (_have_lower) {
1053
        for (int i = 0; i != _arc_num; ++i) {
1054
          Value c = _lower[i];
1055
          if (c >= 0) {
1056
            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
1057
          } else {
1058
            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
1059
          }
1060
          _supply[_source[i]] -= c;
1061
          _supply[_target[i]] += c;
1062
        }
1063
      } else {
1064
        for (int i = 0; i != _arc_num; ++i) {
1065
          _cap[i] = _upper[i];
1066
        }
1074 1067
      }
1075 1068

	
1076
      // Initialize vectors
1077
      _node_num = countNodes(_graph);
1078
      _arc_num = countArcs(_graph);
1079
      int all_node_num = _node_num + 1;
1080
      int all_arc_num = _arc_num + _node_num;
1081
      if (_node_num == 0) return false;
1082

	
1083
      _arc_ref.resize(_arc_num);
1084
      _source.resize(all_arc_num);
1085
      _target.resize(all_arc_num);
1086

	
1087
      _cap.resize(all_arc_num);
1088
      _cost.resize(all_arc_num);
1089
      _supply.resize(all_node_num);
1090
      _flow.resize(all_arc_num);
1091
      _pi.resize(all_node_num);
1092

	
1093
      _parent.resize(all_node_num);
1094
      _pred.resize(all_node_num);
1095
      _forward.resize(all_node_num);
1096
      _thread.resize(all_node_num);
1097
      _rev_thread.resize(all_node_num);
1098
      _succ_num.resize(all_node_num);
1099
      _last_succ.resize(all_node_num);
1100
      _state.resize(all_arc_num);
1101

	
1102
      // Initialize node related data
1103
      bool valid_supply = true;
1104
      Flow sum_supply = 0;
1105
      if (!_pstsup && !_psupply) {
1106
        _pstsup = true;
1107
        _psource = _ptarget = NodeIt(_graph);
1108
        _pstflow = 0;
1109
      }
1110
      if (_psupply) {
1111
        int i = 0;
1112
        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1113
          _node_id[n] = i;
1114
          _supply[i] = (*_psupply)[n];
1115
          sum_supply += _supply[i];
1069
      // Initialize artifical cost
1070
      Cost ART_COST;
1071
      if (std::numeric_limits<Cost>::is_exact) {
1072
        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
1073
      } else {
1074
        ART_COST = std::numeric_limits<Cost>::min();
1075
        for (int i = 0; i != _arc_num; ++i) {
1076
          if (_cost[i] > ART_COST) ART_COST = _cost[i];
1116 1077
        }
1117
        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
1118
                       (_ptype == LEQ && sum_supply >= 0);
1119
      } else {
1120
        int i = 0;
1121
        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1122
          _node_id[n] = i;
1123
          _supply[i] = 0;
1124
        }
1125
        _supply[_node_id[_psource]] =  _pstflow;
1126
        _supply[_node_id[_ptarget]] = -_pstflow;
1127
      }
1128
      if (!valid_supply) return false;
1129

	
1130
      // Infinite capacity value
1131
      Flow inf_cap =
1132
        std::numeric_limits<Flow>::has_infinity ?
1133
        std::numeric_limits<Flow>::infinity() :
1134
        std::numeric_limits<Flow>::max();
1135

	
1136
      // Initialize artifical cost
1137
      Cost art_cost;
1138
      if (std::numeric_limits<Cost>::is_exact) {
1139
        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
1140
      } else {
1141
        art_cost = std::numeric_limits<Cost>::min();
1142
        for (int i = 0; i != _arc_num; ++i) {
1143
          if (_cost[i] > art_cost) art_cost = _cost[i];
1144
        }
1145
        art_cost = (art_cost + 1) * _node_num;
1078
        ART_COST = (ART_COST + 1) * _node_num;
1146 1079
      }
1147 1080

	
1148
      // Run Circulation to check if a feasible solution exists
1149
      typedef ConstMap<Arc, Flow> ConstArcMap;
1150
      ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
1151
      FlowNodeMap *csup = NULL;
1152
      bool local_csup = false;
1153
      if (_psupply) {
1154
        csup = _psupply;
1155
      } else {
1156
        csup = new FlowNodeMap(_graph, 0);
1157
        (*csup)[_psource] =  _pstflow;
1158
        (*csup)[_ptarget] = -_pstflow;
1159
        local_csup = true;
1081
      // Initialize arc maps
1082
      for (int i = 0; i != _arc_num; ++i) {
1083
        _flow[i] = 0;
1084
        _state[i] = STATE_LOWER;
1160 1085
      }
1161
      bool circ_result = false;
1162
      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
1163
        // GEQ problem type
1164
        if (_plower) {
1165
          if (_pupper) {
1166
            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
1167
              circ(_graph, *_plower, *_pupper, *csup);
1168
            circ_result = circ.run();
1169
          } else {
1170
            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
1171
              circ(_graph, *_plower, inf_arc_map, *csup);
1172
            circ_result = circ.run();
1173
          }
1174
        } else {
1175
          if (_pupper) {
1176
            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
1177
              circ(_graph, zero_arc_map, *_pupper, *csup);
1178
            circ_result = circ.run();
1179
          } else {
1180
            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
1181
              circ(_graph, zero_arc_map, inf_arc_map, *csup);
1182
            circ_result = circ.run();
1183
          }
1184
        }
1185
      } else {
1186
        // LEQ problem type
1187
        typedef ReverseDigraph<const GR> RevGraph;
1188
        typedef NegMap<FlowNodeMap> NegNodeMap;
1189
        RevGraph rgraph(_graph);
1190
        NegNodeMap neg_csup(*csup);
1191
        if (_plower) {
1192
          if (_pupper) {
1193
            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
1194
              circ(rgraph, *_plower, *_pupper, neg_csup);
1195
            circ_result = circ.run();
1196
          } else {
1197
            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
1198
              circ(rgraph, *_plower, inf_arc_map, neg_csup);
1199
            circ_result = circ.run();
1200
          }
1201
        } else {
1202
          if (_pupper) {
1203
            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
1204
              circ(rgraph, zero_arc_map, *_pupper, neg_csup);
1205
            circ_result = circ.run();
1206
          } else {
1207
            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
1208
              circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
1209
            circ_result = circ.run();
1210
          }
1211
        }
1212
      }
1213
      if (local_csup) delete csup;
1214
      if (!circ_result) return false;
1215

	
1086
      
1216 1087
      // Set data for the artificial root node
1217 1088
      _root = _node_num;
1218 1089
      _parent[_root] = -1;
1219 1090
      _pred[_root] = -1;
1220 1091
      _thread[_root] = 0;
1221 1092
      _rev_thread[0] = _root;
1222
      _succ_num[_root] = all_node_num;
1093
      _succ_num[_root] = _node_num + 1;
1223 1094
      _last_succ[_root] = _root - 1;
1224
      _supply[_root] = -sum_supply;
1225
      if (sum_supply < 0) {
1226
        _pi[_root] = -art_cost;
1227
      } else {
1228
        _pi[_root] = art_cost;
1229
      }
1230

	
1231
      // Store the arcs in a mixed order
1232
      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1233
      int i = 0;
1234
      for (ArcIt e(_graph); e != INVALID; ++e) {
1235
        _arc_ref[i] = e;
1236
        if ((i += k) >= _arc_num) i = (i % k) + 1;
1237
      }
1238

	
1239
      // Initialize arc maps
1240
      if (_pupper && _pcost) {
1241
        for (int i = 0; i != _arc_num; ++i) {
1242
          Arc e = _arc_ref[i];
1243
          _source[i] = _node_id[_graph.source(e)];
1244
          _target[i] = _node_id[_graph.target(e)];
1245
          _cap[i] = (*_pupper)[e];
1246
          _cost[i] = (*_pcost)[e];
1247
          _flow[i] = 0;
1248
          _state[i] = STATE_LOWER;
1249
        }
1250
      } else {
1251
        for (int i = 0; i != _arc_num; ++i) {
1252
          Arc e = _arc_ref[i];
1253
          _source[i] = _node_id[_graph.source(e)];
1254
          _target[i] = _node_id[_graph.target(e)];
1255
          _flow[i] = 0;
1256
          _state[i] = STATE_LOWER;
1257
        }
1258
        if (_pupper) {
1259
          for (int i = 0; i != _arc_num; ++i)
1260
            _cap[i] = (*_pupper)[_arc_ref[i]];
1261
        } else {
1262
          for (int i = 0; i != _arc_num; ++i)
1263
            _cap[i] = inf_cap;
1264
        }
1265
        if (_pcost) {
1266
          for (int i = 0; i != _arc_num; ++i)
1267
            _cost[i] = (*_pcost)[_arc_ref[i]];
1268
        } else {
1269
          for (int i = 0; i != _arc_num; ++i)
1270
            _cost[i] = 1;
1271
        }
1272
      }
1273
      
1274
      // Remove non-zero lower bounds
1275
      if (_plower) {
1276
        for (int i = 0; i != _arc_num; ++i) {
1277
          Flow c = (*_plower)[_arc_ref[i]];
1278
          if (c != 0) {
1279
            _cap[i] -= c;
1280
            _supply[_source[i]] -= c;
1281
            _supply[_target[i]] += c;
1282
          }
1283
        }
1284
      }
1095
      _supply[_root] = -_sum_supply;
1096
      _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
1285 1097

	
1286 1098
      // Add artificial arcs and initialize the spanning tree data structure
1287 1099
      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1100
        _parent[u] = _root;
1101
        _pred[u] = e;
1288 1102
        _thread[u] = u + 1;
1289 1103
        _rev_thread[u + 1] = u;
1290 1104
        _succ_num[u] = 1;
1291 1105
        _last_succ[u] = u;
1292
        _parent[u] = _root;
1293
        _pred[u] = e;
1294
        _cost[e] = art_cost;
1295
        _cap[e] = inf_cap;
1106
        _cost[e] = ART_COST;
1107
        _cap[e] = INF;
1296 1108
        _state[e] = STATE_TREE;
1297
        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
1109
        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
1298 1110
          _flow[e] = _supply[u];
1299 1111
          _forward[u] = true;
1300
          _pi[u] = -art_cost + _pi[_root];
1112
          _pi[u] = -ART_COST + _pi[_root];
1301 1113
        } else {
1302 1114
          _flow[e] = -_supply[u];
1303 1115
          _forward[u] = false;
1304
          _pi[u] = art_cost + _pi[_root];
1116
          _pi[u] = ART_COST + _pi[_root];
1305 1117
        }
1306 1118
      }
1307 1119

	
1308 1120
      return true;
1309 1121
    }
1310 1122

	
1311 1123
    // Find the join node
1312 1124
    void findJoinNode() {
1313 1125
      int u = _source[in_arc];
1314 1126
      int v = _target[in_arc];
1315 1127
      while (u != v) {
1316 1128
        if (_succ_num[u] < _succ_num[v]) {
1317 1129
          u = _parent[u];
1318 1130
        } else {
1319 1131
          v = _parent[v];
1320 1132
        }
1321 1133
      }
1322 1134
      join = u;
1323 1135
    }
1324 1136

	
1325 1137
    // Find the leaving arc of the cycle and returns true if the
1326 1138
    // leaving arc is not the same as the entering arc
1327 1139
    bool findLeavingArc() {
1328 1140
      // Initialize first and second nodes according to the direction
1329 1141
      // of the cycle
1330 1142
      if (_state[in_arc] == STATE_LOWER) {
1331 1143
        first  = _source[in_arc];
1332 1144
        second = _target[in_arc];
1333 1145
      } else {
1334 1146
        first  = _target[in_arc];
1335 1147
        second = _source[in_arc];
1336 1148
      }
1337 1149
      delta = _cap[in_arc];
1338 1150
      int result = 0;
1339
      Flow d;
1151
      Value d;
1340 1152
      int e;
1341 1153

	
1342 1154
      // Search the cycle along the path form the first node to the root
1343 1155
      for (int u = first; u != join; u = _parent[u]) {
1344 1156
        e = _pred[u];
1345
        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
1157
        d = _forward[u] ?
1158
          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
1346 1159
        if (d < delta) {
1347 1160
          delta = d;
1348 1161
          u_out = u;
1349 1162
          result = 1;
1350 1163
        }
1351 1164
      }
1352 1165
      // Search the cycle along the path form the second node to the root
1353 1166
      for (int u = second; u != join; u = _parent[u]) {
1354 1167
        e = _pred[u];
1355
        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
1168
        d = _forward[u] ? 
1169
          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
1356 1170
        if (d <= delta) {
1357 1171
          delta = d;
1358 1172
          u_out = u;
1359 1173
          result = 2;
1360 1174
        }
1361 1175
      }
1362 1176

	
1363 1177
      if (result == 1) {
1364 1178
        u_in = first;
1365 1179
        v_in = second;
1366 1180
      } else {
1367 1181
        u_in = second;
1368 1182
        v_in = first;
1369 1183
      }
1370 1184
      return result != 0;
1371 1185
    }
1372 1186

	
1373 1187
    // Change _flow and _state vectors
1374 1188
    void changeFlow(bool change) {
1375 1189
      // Augment along the cycle
1376 1190
      if (delta > 0) {
1377
        Flow val = _state[in_arc] * delta;
1191
        Value val = _state[in_arc] * delta;
1378 1192
        _flow[in_arc] += val;
1379 1193
        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
1380 1194
          _flow[_pred[u]] += _forward[u] ? -val : val;
1381 1195
        }
1382 1196
        for (int u = _target[in_arc]; u != join; u = _parent[u]) {
1383 1197
          _flow[_pred[u]] += _forward[u] ? val : -val;
1384 1198
        }
1385 1199
      }
1386 1200
      // Update the state of the entering and leaving arcs
1387 1201
      if (change) {
1388 1202
        _state[in_arc] = STATE_TREE;
1389 1203
        _state[_pred[u_out]] =
1390 1204
          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
1391 1205
      } else {
1392 1206
        _state[in_arc] = -_state[in_arc];
1393 1207
      }
1394 1208
    }
1395 1209

	
1396 1210
    // Update the tree structure
1397 1211
    void updateTreeStructure() {
1398 1212
      int u, w;
1399 1213
      int old_rev_thread = _rev_thread[u_out];
1400 1214
      int old_succ_num = _succ_num[u_out];
1401 1215
      int old_last_succ = _last_succ[u_out];
1402 1216
      v_out = _parent[u_out];
1403 1217

	
1404 1218
      u = _last_succ[u_in];  // the last successor of u_in
1405 1219
      right = _thread[u];    // the node after it
1406 1220

	
1407 1221
      // Handle the case when old_rev_thread equals to v_in
1408 1222
      // (it also means that join and v_out coincide)
1409 1223
      if (old_rev_thread == v_in) {
1410 1224
        last = _thread[_last_succ[u_out]];
1411 1225
      } else {
1412 1226
        last = _thread[v_in];
1413 1227
      }
1414 1228

	
1415 1229
      // Update _thread and _parent along the stem nodes (i.e. the nodes
1416 1230
      // between u_in and u_out, whose parent have to be changed)
1417 1231
      _thread[v_in] = stem = u_in;
1418 1232
      _dirty_revs.clear();
1419 1233
      _dirty_revs.push_back(v_in);
1420 1234
      par_stem = v_in;
1421 1235
      while (stem != u_out) {
1422 1236
        // Insert the next stem node into the thread list
1423 1237
        new_stem = _parent[stem];
1424 1238
        _thread[u] = new_stem;
1425 1239
        _dirty_revs.push_back(u);
1426 1240

	
1427 1241
        // Remove the subtree of stem from the thread list
1428 1242
        w = _rev_thread[stem];
1429 1243
        _thread[w] = right;
1430 1244
        _rev_thread[right] = w;
1431 1245

	
1432 1246
        // Change the parent node and shift stem nodes
1433 1247
        _parent[stem] = par_stem;
1434 1248
        par_stem = stem;
1435 1249
        stem = new_stem;
1436 1250

	
1437 1251
        // Update u and right
1438 1252
        u = _last_succ[stem] == _last_succ[par_stem] ?
1439 1253
          _rev_thread[par_stem] : _last_succ[stem];
1440 1254
        right = _thread[u];
1441 1255
      }
1442 1256
      _parent[u_out] = par_stem;
1443 1257
      _thread[u] = last;
1444 1258
      _rev_thread[last] = u;
1445 1259
      _last_succ[u_out] = u;
1446 1260

	
1447 1261
      // Remove the subtree of u_out from the thread list except for
1448 1262
      // the case when old_rev_thread equals to v_in
1449 1263
      // (it also means that join and v_out coincide)
1450 1264
      if (old_rev_thread != v_in) {
1451 1265
        _thread[old_rev_thread] = right;
1452 1266
        _rev_thread[right] = old_rev_thread;
1453 1267
      }
1454 1268

	
1455 1269
      // Update _rev_thread using the new _thread values
1456 1270
      for (int i = 0; i < int(_dirty_revs.size()); ++i) {
1457 1271
        u = _dirty_revs[i];
1458 1272
        _rev_thread[_thread[u]] = u;
1459 1273
      }
1460 1274

	
1461 1275
      // Update _pred, _forward, _last_succ and _succ_num for the
1462 1276
      // stem nodes from u_out to u_in
1463 1277
      int tmp_sc = 0, tmp_ls = _last_succ[u_out];
1464 1278
      u = u_out;
1465 1279
      while (u != u_in) {
1466 1280
        w = _parent[u];
1467 1281
        _pred[u] = _pred[w];
1468 1282
        _forward[u] = !_forward[w];
1469 1283
        tmp_sc += _succ_num[u] - _succ_num[w];
1470 1284
        _succ_num[u] = tmp_sc;
1471 1285
        _last_succ[w] = tmp_ls;
1472 1286
        u = w;
1473 1287
      }
1474 1288
      _pred[u_in] = in_arc;
1475 1289
      _forward[u_in] = (u_in == _source[in_arc]);
1476 1290
      _succ_num[u_in] = old_succ_num;
1477 1291

	
1478 1292
      // Set limits for updating _last_succ form v_in and v_out
1479 1293
      // towards the root
1480 1294
      int up_limit_in = -1;
1481 1295
      int up_limit_out = -1;
1482 1296
      if (_last_succ[join] == v_in) {
1483 1297
        up_limit_out = join;
1484 1298
      } else {
1485 1299
        up_limit_in = join;
1486 1300
      }
1487 1301

	
1488 1302
      // Update _last_succ from v_in towards the root
1489 1303
      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
1490 1304
           u = _parent[u]) {
1491 1305
        _last_succ[u] = _last_succ[u_out];
1492 1306
      }
1493 1307
      // Update _last_succ from v_out towards the root
1494 1308
      if (join != old_rev_thread && v_in != old_rev_thread) {
1495 1309
        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1496 1310
             u = _parent[u]) {
1497 1311
          _last_succ[u] = old_rev_thread;
1498 1312
        }
1499 1313
      } else {
1500 1314
        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1501 1315
             u = _parent[u]) {
1502 1316
          _last_succ[u] = _last_succ[u_out];
1503 1317
        }
1504 1318
      }
1505 1319

	
1506 1320
      // Update _succ_num from v_in to join
1507 1321
      for (u = v_in; u != join; u = _parent[u]) {
1508 1322
        _succ_num[u] += old_succ_num;
1509 1323
      }
1510 1324
      // Update _succ_num from v_out to join
1511 1325
      for (u = v_out; u != join; u = _parent[u]) {
1512 1326
        _succ_num[u] -= old_succ_num;
1513 1327
      }
1514 1328
    }
1515 1329

	
1516 1330
    // Update potentials
1517 1331
    void updatePotential() {
1518 1332
      Cost sigma = _forward[u_in] ?
1519 1333
        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1520 1334
        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1521 1335
      // Update potentials in the subtree, which has been moved
1522 1336
      int end = _thread[_last_succ[u_in]];
1523 1337
      for (int u = u_in; u != end; u = _thread[u]) {
1524 1338
        _pi[u] += sigma;
1525 1339
      }
1526 1340
    }
1527 1341

	
1528 1342
    // Execute the algorithm
1529
    bool start(PivotRule pivot_rule) {
1343
    ProblemType start(PivotRule pivot_rule) {
1530 1344
      // Select the pivot rule implementation
1531 1345
      switch (pivot_rule) {
1532 1346
        case FIRST_ELIGIBLE:
1533 1347
          return start<FirstEligiblePivotRule>();
1534 1348
        case BEST_ELIGIBLE:
1535 1349
          return start<BestEligiblePivotRule>();
1536 1350
        case BLOCK_SEARCH:
1537 1351
          return start<BlockSearchPivotRule>();
1538 1352
        case CANDIDATE_LIST:
1539 1353
          return start<CandidateListPivotRule>();
1540 1354
        case ALTERING_LIST:
1541 1355
          return start<AlteringListPivotRule>();
1542 1356
      }
1543
      return false;
1357
      return INFEASIBLE; // avoid warning
1544 1358
    }
1545 1359

	
1546 1360
    template <typename PivotRuleImpl>
1547
    bool start() {
1361
    ProblemType start() {
1548 1362
      PivotRuleImpl pivot(*this);
1549 1363

	
1550 1364
      // Execute the Network Simplex algorithm
1551 1365
      while (pivot.findEnteringArc()) {
1552 1366
        findJoinNode();
1553 1367
        bool change = findLeavingArc();
1368
        if (delta >= INF) return UNBOUNDED;
1554 1369
        changeFlow(change);
1555 1370
        if (change) {
1556 1371
          updateTreeStructure();
1557 1372
          updatePotential();
1558 1373
        }
1559 1374
      }
1560

	
1561
      // Copy flow values to _flow_map
1562
      if (_plower) {
1563
        for (int i = 0; i != _arc_num; ++i) {
1564
          Arc e = _arc_ref[i];
1565
          _flow_map->set(e, (*_plower)[e] + _flow[i]);
1566
        }
1567
      } else {
1568
        for (int i = 0; i != _arc_num; ++i) {
1569
          _flow_map->set(_arc_ref[i], _flow[i]);
1375
      
1376
      // Check feasibility
1377
      if (_sum_supply < 0) {
1378
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1379
          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
1570 1380
        }
1571 1381
      }
1572
      // Copy potential values to _potential_map
1573
      for (NodeIt n(_graph); n != INVALID; ++n) {
1574
        _potential_map->set(n, _pi[_node_id[n]]);
1382
      else if (_sum_supply > 0) {
1383
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1384
          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
1385
        }
1386
      }
1387
      else {
1388
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1389
          if (_flow[e] != 0) return INFEASIBLE;
1390
        }
1575 1391
      }
1576 1392

	
1577
      return true;
1393
      // Transform the solution and the supply map to the original form
1394
      if (_have_lower) {
1395
        for (int i = 0; i != _arc_num; ++i) {
1396
          Value c = _lower[i];
1397
          if (c != 0) {
1398
            _flow[i] += c;
1399
            _supply[_source[i]] += c;
1400
            _supply[_target[i]] -= c;
1401
          }
1402
        }
1403
      }
1404

	
1405
      return OPTIMAL;
1578 1406
    }
1579 1407

	
1580 1408
  }; //class NetworkSimplex
1581 1409

	
1582 1410
  ///@}
1583 1411

	
1584 1412
} //namespace lemon
1585 1413

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

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

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

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

	
29 29
namespace lemon {
30 30

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

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

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

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

	
51 51
    /// \brief The type of the map that stores the flow values.
52 52
    ///
53 53
    /// The type of the map that stores the flow values.
54 54
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
55
    typedef typename Digraph::template ArcMap<Flow> FlowMap;
55
    typedef typename Digraph::template ArcMap<Value> FlowMap;
56 56

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

	
66 66
    /// \brief The elevator type used by Preflow algorithm.
67 67
    ///
68 68
    /// The elevator type used by Preflow algorithm.
69 69
    ///
70 70
    /// \sa Elevator
71 71
    /// \sa LinkedElevator
72 72
    typedef LinkedElevator<Digraph, typename Digraph::Node> Elevator;
73 73

	
74 74
    /// \brief Instantiates an Elevator.
75 75
    ///
76 76
    /// This function instantiates an \ref Elevator.
77 77
    /// \param digraph The digraph for which we would like to define
78 78
    /// the elevator.
79 79
    /// \param max_level The maximum level of the elevator.
80 80
    static Elevator* createElevator(const Digraph& digraph, int max_level) {
81 81
      return new Elevator(digraph, max_level);
82 82
    }
83 83

	
84 84
    /// \brief The tolerance used by the algorithm
85 85
    ///
86 86
    /// The tolerance used by the algorithm to handle inexact computation.
87
    typedef lemon::Tolerance<Flow> Tolerance;
87
    typedef lemon::Tolerance<Value> Tolerance;
88 88

	
89 89
  };
90 90

	
91 91

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

	
121 121
    ///The \ref PreflowDefaultTraits "traits class" of the algorithm.
122 122
    typedef TR Traits;
123 123
    ///The type of the digraph the algorithm runs on.
124 124
    typedef typename Traits::Digraph Digraph;
125 125
    ///The type of the capacity map.
126 126
    typedef typename Traits::CapacityMap CapacityMap;
127 127
    ///The type of the flow values.
128
    typedef typename Traits::Flow Flow;
128
    typedef typename Traits::Value Value;
129 129

	
130 130
    ///The type of the flow map.
131 131
    typedef typename Traits::FlowMap FlowMap;
132 132
    ///The type of the elevator.
133 133
    typedef typename Traits::Elevator Elevator;
134 134
    ///The type of the tolerance.
135 135
    typedef typename Traits::Tolerance Tolerance;
136 136

	
137 137
  private:
138 138

	
139 139
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
140 140

	
141 141
    const Digraph& _graph;
142 142
    const CapacityMap* _capacity;
143 143

	
144 144
    int _node_num;
145 145

	
146 146
    Node _source, _target;
147 147

	
148 148
    FlowMap* _flow;
149 149
    bool _local_flow;
150 150

	
151 151
    Elevator* _level;
152 152
    bool _local_level;
153 153

	
154
    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
154
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
155 155
    ExcessMap* _excess;
156 156

	
157 157
    Tolerance _tolerance;
158 158

	
159 159
    bool _phase;
160 160

	
161 161

	
162 162
    void createStructures() {
163 163
      _node_num = countNodes(_graph);
164 164

	
165 165
      if (!_flow) {
166 166
        _flow = Traits::createFlowMap(_graph);
167 167
        _local_flow = true;
168 168
      }
169 169
      if (!_level) {
170 170
        _level = Traits::createElevator(_graph, _node_num);
171 171
        _local_level = true;
172 172
      }
173 173
      if (!_excess) {
174 174
        _excess = new ExcessMap(_graph);
175 175
      }
176 176
    }
177 177

	
178 178
    void destroyStructures() {
179 179
      if (_local_flow) {
180 180
        delete _flow;
181 181
      }
182 182
      if (_local_level) {
183 183
        delete _level;
184 184
      }
185 185
      if (_excess) {
186 186
        delete _excess;
187 187
      }
188 188
    }
189 189

	
190 190
  public:
191 191

	
192 192
    typedef Preflow Create;
193 193

	
194 194
    ///\name Named Template Parameters
195 195

	
196 196
    ///@{
197 197

	
198 198
    template <typename T>
199 199
    struct SetFlowMapTraits : public Traits {
200 200
      typedef T FlowMap;
201 201
      static FlowMap *createFlowMap(const Digraph&) {
202 202
        LEMON_ASSERT(false, "FlowMap is not initialized");
203 203
        return 0; // ignore warnings
204 204
      }
205 205
    };
206 206

	
207 207
    /// \brief \ref named-templ-param "Named parameter" for setting
208 208
    /// FlowMap type
209 209
    ///
210 210
    /// \ref named-templ-param "Named parameter" for setting FlowMap
211 211
    /// type.
212 212
    template <typename T>
213 213
    struct SetFlowMap
214 214
      : public Preflow<Digraph, CapacityMap, SetFlowMapTraits<T> > {
215 215
      typedef Preflow<Digraph, CapacityMap,
216 216
                      SetFlowMapTraits<T> > Create;
217 217
    };
218 218

	
219 219
    template <typename T>
220 220
    struct SetElevatorTraits : public Traits {
221 221
      typedef T Elevator;
222 222
      static Elevator *createElevator(const Digraph&, int) {
223 223
        LEMON_ASSERT(false, "Elevator is not initialized");
224 224
        return 0; // ignore warnings
225 225
      }
226 226
    };
227 227

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

	
244 244
    template <typename T>
245 245
    struct SetStandardElevatorTraits : public Traits {
246 246
      typedef T Elevator;
247 247
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
248 248
        return new Elevator(digraph, max_level);
249 249
      }
250 250
    };
251 251

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

	
272 272
    /// @}
273 273

	
274 274
  protected:
275 275

	
276 276
    Preflow() {}
277 277

	
278 278
  public:
279 279

	
280 280

	
281 281
    /// \brief The constructor of the class.
282 282
    ///
283 283
    /// The constructor of the class.
284 284
    /// \param digraph The digraph the algorithm runs on.
285 285
    /// \param capacity The capacity of the arcs.
286 286
    /// \param source The source node.
287 287
    /// \param target The target node.
288 288
    Preflow(const Digraph& digraph, const CapacityMap& capacity,
289 289
            Node source, Node target)
290 290
      : _graph(digraph), _capacity(&capacity),
291 291
        _node_num(0), _source(source), _target(target),
292 292
        _flow(0), _local_flow(false),
293 293
        _level(0), _local_level(false),
294 294
        _excess(0), _tolerance(), _phase() {}
295 295

	
296 296
    /// \brief Destructor.
297 297
    ///
298 298
    /// Destructor.
299 299
    ~Preflow() {
300 300
      destroyStructures();
301 301
    }
302 302

	
303 303
    /// \brief Sets the capacity map.
304 304
    ///
305 305
    /// Sets the capacity map.
306 306
    /// \return <tt>(*this)</tt>
307 307
    Preflow& capacityMap(const CapacityMap& map) {
308 308
      _capacity = &map;
309 309
      return *this;
310 310
    }
311 311

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

	
329 329
    /// \brief Sets the source node.
330 330
    ///
331 331
    /// Sets the source node.
332 332
    /// \return <tt>(*this)</tt>
333 333
    Preflow& source(const Node& node) {
334 334
      _source = node;
335 335
      return *this;
336 336
    }
337 337

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

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

	
364 364
    /// \brief Returns a const reference to the elevator.
365 365
    ///
366 366
    /// Returns a const reference to the elevator.
367 367
    ///
368 368
    /// \pre Either \ref run() or \ref init() must be called before
369 369
    /// using this function.
370 370
    const Elevator& elevator() const {
371 371
      return *_level;
372 372
    }
373 373

	
374 374
    /// \brief Sets the tolerance used by algorithm.
375 375
    ///
376 376
    /// Sets the tolerance used by algorithm.
377 377
    Preflow& tolerance(const Tolerance& tolerance) const {
378 378
      _tolerance = tolerance;
379 379
      return *this;
380 380
    }
381 381

	
382 382
    /// \brief Returns a const reference to the tolerance.
383 383
    ///
384 384
    /// Returns a const reference to the tolerance.
385 385
    const Tolerance& tolerance() const {
386 386
      return tolerance;
387 387
    }
388 388

	
389 389
    /// \name Execution Control
390 390
    /// The simplest way to execute the preflow algorithm is to use
391 391
    /// \ref run() or \ref runMinCut().\n
392 392
    /// If you need more control on the initial solution or the execution,
393 393
    /// first you have to call one of the \ref init() functions, then
394 394
    /// \ref startFirstPhase() and if you need it \ref startSecondPhase().
395 395

	
396 396
    ///@{
397 397

	
398 398
    /// \brief Initializes the internal data structures.
399 399
    ///
400 400
    /// Initializes the internal data structures and sets the initial
401 401
    /// flow to zero on each arc.
402 402
    void init() {
403 403
      createStructures();
404 404

	
405 405
      _phase = true;
406 406
      for (NodeIt n(_graph); n != INVALID; ++n) {
407 407
        (*_excess)[n] = 0;
408 408
      }
409 409

	
410 410
      for (ArcIt e(_graph); e != INVALID; ++e) {
411 411
        _flow->set(e, 0);
412 412
      }
413 413

	
414 414
      typename Digraph::template NodeMap<bool> reached(_graph, false);
415 415

	
416 416
      _level->initStart();
417 417
      _level->initAddItem(_target);
418 418

	
419 419
      std::vector<Node> queue;
420 420
      reached[_source] = true;
421 421

	
422 422
      queue.push_back(_target);
423 423
      reached[_target] = true;
424 424
      while (!queue.empty()) {
425 425
        _level->initNewLevel();
426 426
        std::vector<Node> nqueue;
427 427
        for (int i = 0; i < int(queue.size()); ++i) {
428 428
          Node n = queue[i];
429 429
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
430 430
            Node u = _graph.source(e);
431 431
            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
432 432
              reached[u] = true;
433 433
              _level->initAddItem(u);
434 434
              nqueue.push_back(u);
435 435
            }
436 436
          }
437 437
        }
438 438
        queue.swap(nqueue);
439 439
      }
440 440
      _level->initFinish();
441 441

	
442 442
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
443 443
        if (_tolerance.positive((*_capacity)[e])) {
444 444
          Node u = _graph.target(e);
445 445
          if ((*_level)[u] == _level->maxLevel()) continue;
446 446
          _flow->set(e, (*_capacity)[e]);
447 447
          (*_excess)[u] += (*_capacity)[e];
448 448
          if (u != _target && !_level->active(u)) {
449 449
            _level->activate(u);
450 450
          }
451 451
        }
452 452
      }
453 453
    }
454 454

	
455 455
    /// \brief Initializes the internal data structures using the
456 456
    /// given flow map.
457 457
    ///
458 458
    /// Initializes the internal data structures and sets the initial
459 459
    /// flow to the given \c flowMap. The \c flowMap should contain a
460 460
    /// flow or at least a preflow, i.e. at each node excluding the
461 461
    /// source node the incoming flow should greater or equal to the
462 462
    /// outgoing flow.
463 463
    /// \return \c false if the given \c flowMap is not a preflow.
464 464
    template <typename FlowMap>
465 465
    bool init(const FlowMap& flowMap) {
466 466
      createStructures();
467 467

	
468 468
      for (ArcIt e(_graph); e != INVALID; ++e) {
469 469
        _flow->set(e, flowMap[e]);
470 470
      }
471 471

	
472 472
      for (NodeIt n(_graph); n != INVALID; ++n) {
473
        Flow excess = 0;
473
        Value excess = 0;
474 474
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
475 475
          excess += (*_flow)[e];
476 476
        }
477 477
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
478 478
          excess -= (*_flow)[e];
479 479
        }
480 480
        if (excess < 0 && n != _source) return false;
481 481
        (*_excess)[n] = excess;
482 482
      }
483 483

	
484 484
      typename Digraph::template NodeMap<bool> reached(_graph, false);
485 485

	
486 486
      _level->initStart();
487 487
      _level->initAddItem(_target);
488 488

	
489 489
      std::vector<Node> queue;
490 490
      reached[_source] = true;
491 491

	
492 492
      queue.push_back(_target);
493 493
      reached[_target] = true;
494 494
      while (!queue.empty()) {
495 495
        _level->initNewLevel();
496 496
        std::vector<Node> nqueue;
497 497
        for (int i = 0; i < int(queue.size()); ++i) {
498 498
          Node n = queue[i];
499 499
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
500 500
            Node u = _graph.source(e);
501 501
            if (!reached[u] &&
502 502
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
503 503
              reached[u] = true;
504 504
              _level->initAddItem(u);
505 505
              nqueue.push_back(u);
506 506
            }
507 507
          }
508 508
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
509 509
            Node v = _graph.target(e);
510 510
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
511 511
              reached[v] = true;
512 512
              _level->initAddItem(v);
513 513
              nqueue.push_back(v);
514 514
            }
515 515
          }
516 516
        }
517 517
        queue.swap(nqueue);
518 518
      }
519 519
      _level->initFinish();
520 520

	
521 521
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
522
        Flow rem = (*_capacity)[e] - (*_flow)[e];
522
        Value rem = (*_capacity)[e] - (*_flow)[e];
523 523
        if (_tolerance.positive(rem)) {
524 524
          Node u = _graph.target(e);
525 525
          if ((*_level)[u] == _level->maxLevel()) continue;
526 526
          _flow->set(e, (*_capacity)[e]);
527 527
          (*_excess)[u] += rem;
528 528
          if (u != _target && !_level->active(u)) {
529 529
            _level->activate(u);
530 530
          }
531 531
        }
532 532
      }
533 533
      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
534
        Flow rem = (*_flow)[e];
534
        Value rem = (*_flow)[e];
535 535
        if (_tolerance.positive(rem)) {
536 536
          Node v = _graph.source(e);
537 537
          if ((*_level)[v] == _level->maxLevel()) continue;
538 538
          _flow->set(e, 0);
539 539
          (*_excess)[v] += rem;
540 540
          if (v != _target && !_level->active(v)) {
541 541
            _level->activate(v);
542 542
          }
543 543
        }
544 544
      }
545 545
      return true;
546 546
    }
547 547

	
548 548
    /// \brief Starts the first phase of the preflow algorithm.
549 549
    ///
550 550
    /// The preflow algorithm consists of two phases, this method runs
551 551
    /// the first phase. After the first phase the maximum flow value
552 552
    /// and a minimum value cut can already be computed, although a
553 553
    /// maximum flow is not yet obtained. So after calling this method
554 554
    /// \ref flowValue() returns the value of a maximum flow and \ref
555 555
    /// minCut() returns a minimum cut.
556 556
    /// \pre One of the \ref init() functions must be called before
557 557
    /// using this function.
558 558
    void startFirstPhase() {
559 559
      _phase = true;
560 560

	
561 561
      Node n = _level->highestActive();
562 562
      int level = _level->highestActiveLevel();
563 563
      while (n != INVALID) {
564 564
        int num = _node_num;
565 565

	
566 566
        while (num > 0 && n != INVALID) {
567
          Flow excess = (*_excess)[n];
567
          Value excess = (*_excess)[n];
568 568
          int new_level = _level->maxLevel();
569 569

	
570 570
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
571
            Flow rem = (*_capacity)[e] - (*_flow)[e];
571
            Value rem = (*_capacity)[e] - (*_flow)[e];
572 572
            if (!_tolerance.positive(rem)) continue;
573 573
            Node v = _graph.target(e);
574 574
            if ((*_level)[v] < level) {
575 575
              if (!_level->active(v) && v != _target) {
576 576
                _level->activate(v);
577 577
              }
578 578
              if (!_tolerance.less(rem, excess)) {
579 579
                _flow->set(e, (*_flow)[e] + excess);
580 580
                (*_excess)[v] += excess;
581 581
                excess = 0;
582 582
                goto no_more_push_1;
583 583
              } else {
584 584
                excess -= rem;
585 585
                (*_excess)[v] += rem;
586 586
                _flow->set(e, (*_capacity)[e]);
587 587
              }
588 588
            } else if (new_level > (*_level)[v]) {
589 589
              new_level = (*_level)[v];
590 590
            }
591 591
          }
592 592

	
593 593
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
594
            Flow rem = (*_flow)[e];
594
            Value rem = (*_flow)[e];
595 595
            if (!_tolerance.positive(rem)) continue;
596 596
            Node v = _graph.source(e);
597 597
            if ((*_level)[v] < level) {
598 598
              if (!_level->active(v) && v != _target) {
599 599
                _level->activate(v);
600 600
              }
601 601
              if (!_tolerance.less(rem, excess)) {
602 602
                _flow->set(e, (*_flow)[e] - excess);
603 603
                (*_excess)[v] += excess;
604 604
                excess = 0;
605 605
                goto no_more_push_1;
606 606
              } else {
607 607
                excess -= rem;
608 608
                (*_excess)[v] += rem;
609 609
                _flow->set(e, 0);
610 610
              }
611 611
            } else if (new_level > (*_level)[v]) {
612 612
              new_level = (*_level)[v];
613 613
            }
614 614
          }
615 615

	
616 616
        no_more_push_1:
617 617

	
618 618
          (*_excess)[n] = excess;
619 619

	
620 620
          if (excess != 0) {
621 621
            if (new_level + 1 < _level->maxLevel()) {
622 622
              _level->liftHighestActive(new_level + 1);
623 623
            } else {
624 624
              _level->liftHighestActiveToTop();
625 625
            }
626 626
            if (_level->emptyLevel(level)) {
627 627
              _level->liftToTop(level);
628 628
            }
629 629
          } else {
630 630
            _level->deactivate(n);
631 631
          }
632 632

	
633 633
          n = _level->highestActive();
634 634
          level = _level->highestActiveLevel();
635 635
          --num;
636 636
        }
637 637

	
638 638
        num = _node_num * 20;
639 639
        while (num > 0 && n != INVALID) {
640
          Flow excess = (*_excess)[n];
640
          Value excess = (*_excess)[n];
641 641
          int new_level = _level->maxLevel();
642 642

	
643 643
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
644
            Flow rem = (*_capacity)[e] - (*_flow)[e];
644
            Value rem = (*_capacity)[e] - (*_flow)[e];
645 645
            if (!_tolerance.positive(rem)) continue;
646 646
            Node v = _graph.target(e);
647 647
            if ((*_level)[v] < level) {
648 648
              if (!_level->active(v) && v != _target) {
649 649
                _level->activate(v);
650 650
              }
651 651
              if (!_tolerance.less(rem, excess)) {
652 652
                _flow->set(e, (*_flow)[e] + excess);
653 653
                (*_excess)[v] += excess;
654 654
                excess = 0;
655 655
                goto no_more_push_2;
656 656
              } else {
657 657
                excess -= rem;
658 658
                (*_excess)[v] += rem;
659 659
                _flow->set(e, (*_capacity)[e]);
660 660
              }
661 661
            } else if (new_level > (*_level)[v]) {
662 662
              new_level = (*_level)[v];
663 663
            }
664 664
          }
665 665

	
666 666
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
667
            Flow rem = (*_flow)[e];
667
            Value rem = (*_flow)[e];
668 668
            if (!_tolerance.positive(rem)) continue;
669 669
            Node v = _graph.source(e);
670 670
            if ((*_level)[v] < level) {
671 671
              if (!_level->active(v) && v != _target) {
672 672
                _level->activate(v);
673 673
              }
674 674
              if (!_tolerance.less(rem, excess)) {
675 675
                _flow->set(e, (*_flow)[e] - excess);
676 676
                (*_excess)[v] += excess;
677 677
                excess = 0;
678 678
                goto no_more_push_2;
679 679
              } else {
680 680
                excess -= rem;
681 681
                (*_excess)[v] += rem;
682 682
                _flow->set(e, 0);
683 683
              }
684 684
            } else if (new_level > (*_level)[v]) {
685 685
              new_level = (*_level)[v];
686 686
            }
687 687
          }
688 688

	
689 689
        no_more_push_2:
690 690

	
691 691
          (*_excess)[n] = excess;
692 692

	
693 693
          if (excess != 0) {
694 694
            if (new_level + 1 < _level->maxLevel()) {
695 695
              _level->liftActiveOn(level, new_level + 1);
696 696
            } else {
697 697
              _level->liftActiveToTop(level);
698 698
            }
699 699
            if (_level->emptyLevel(level)) {
700 700
              _level->liftToTop(level);
701 701
            }
702 702
          } else {
703 703
            _level->deactivate(n);
704 704
          }
705 705

	
706 706
          while (level >= 0 && _level->activeFree(level)) {
707 707
            --level;
708 708
          }
709 709
          if (level == -1) {
710 710
            n = _level->highestActive();
711 711
            level = _level->highestActiveLevel();
712 712
          } else {
713 713
            n = _level->activeOn(level);
714 714
          }
715 715
          --num;
716 716
        }
717 717
      }
718 718
    }
719 719

	
720 720
    /// \brief Starts the second phase of the preflow algorithm.
721 721
    ///
722 722
    /// The preflow algorithm consists of two phases, this method runs
723 723
    /// the second phase. After calling one of the \ref init() functions
724 724
    /// and \ref startFirstPhase() and then \ref startSecondPhase(),
725 725
    /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
726 726
    /// value of a maximum flow, \ref minCut() returns a minimum cut
727 727
    /// \pre One of the \ref init() functions and \ref startFirstPhase()
728 728
    /// must be called before using this function.
729 729
    void startSecondPhase() {
730 730
      _phase = false;
731 731

	
732 732
      typename Digraph::template NodeMap<bool> reached(_graph);
733 733
      for (NodeIt n(_graph); n != INVALID; ++n) {
734 734
        reached[n] = (*_level)[n] < _level->maxLevel();
735 735
      }
736 736

	
737 737
      _level->initStart();
738 738
      _level->initAddItem(_source);
739 739

	
740 740
      std::vector<Node> queue;
741 741
      queue.push_back(_source);
742 742
      reached[_source] = true;
743 743

	
744 744
      while (!queue.empty()) {
745 745
        _level->initNewLevel();
746 746
        std::vector<Node> nqueue;
747 747
        for (int i = 0; i < int(queue.size()); ++i) {
748 748
          Node n = queue[i];
749 749
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
750 750
            Node v = _graph.target(e);
751 751
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
752 752
              reached[v] = true;
753 753
              _level->initAddItem(v);
754 754
              nqueue.push_back(v);
755 755
            }
756 756
          }
757 757
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
758 758
            Node u = _graph.source(e);
759 759
            if (!reached[u] &&
760 760
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
761 761
              reached[u] = true;
762 762
              _level->initAddItem(u);
763 763
              nqueue.push_back(u);
764 764
            }
765 765
          }
766 766
        }
767 767
        queue.swap(nqueue);
768 768
      }
769 769
      _level->initFinish();
770 770

	
771 771
      for (NodeIt n(_graph); n != INVALID; ++n) {
772 772
        if (!reached[n]) {
773 773
          _level->dirtyTopButOne(n);
774 774
        } else if ((*_excess)[n] > 0 && _target != n) {
775 775
          _level->activate(n);
776 776
        }
777 777
      }
778 778

	
779 779
      Node n;
780 780
      while ((n = _level->highestActive()) != INVALID) {
781
        Flow excess = (*_excess)[n];
781
        Value excess = (*_excess)[n];
782 782
        int level = _level->highestActiveLevel();
783 783
        int new_level = _level->maxLevel();
784 784

	
785 785
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
786
          Flow rem = (*_capacity)[e] - (*_flow)[e];
786
          Value rem = (*_capacity)[e] - (*_flow)[e];
787 787
          if (!_tolerance.positive(rem)) continue;
788 788
          Node v = _graph.target(e);
789 789
          if ((*_level)[v] < level) {
790 790
            if (!_level->active(v) && v != _source) {
791 791
              _level->activate(v);
792 792
            }
793 793
            if (!_tolerance.less(rem, excess)) {
794 794
              _flow->set(e, (*_flow)[e] + excess);
795 795
              (*_excess)[v] += excess;
796 796
              excess = 0;
797 797
              goto no_more_push;
798 798
            } else {
799 799
              excess -= rem;
800 800
              (*_excess)[v] += rem;
801 801
              _flow->set(e, (*_capacity)[e]);
802 802
            }
803 803
          } else if (new_level > (*_level)[v]) {
804 804
            new_level = (*_level)[v];
805 805
          }
806 806
        }
807 807

	
808 808
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
809
          Flow rem = (*_flow)[e];
809
          Value rem = (*_flow)[e];
810 810
          if (!_tolerance.positive(rem)) continue;
811 811
          Node v = _graph.source(e);
812 812
          if ((*_level)[v] < level) {
813 813
            if (!_level->active(v) && v != _source) {
814 814
              _level->activate(v);
815 815
            }
816 816
            if (!_tolerance.less(rem, excess)) {
817 817
              _flow->set(e, (*_flow)[e] - excess);
818 818
              (*_excess)[v] += excess;
819 819
              excess = 0;
820 820
              goto no_more_push;
821 821
            } else {
822 822
              excess -= rem;
823 823
              (*_excess)[v] += rem;
824 824
              _flow->set(e, 0);
825 825
            }
826 826
          } else if (new_level > (*_level)[v]) {
827 827
            new_level = (*_level)[v];
828 828
          }
829 829
        }
830 830

	
831 831
      no_more_push:
832 832

	
833 833
        (*_excess)[n] = excess;
834 834

	
835 835
        if (excess != 0) {
836 836
          if (new_level + 1 < _level->maxLevel()) {
837 837
            _level->liftHighestActive(new_level + 1);
838 838
          } else {
839 839
            // Calculation error
840 840
            _level->liftHighestActiveToTop();
841 841
          }
842 842
          if (_level->emptyLevel(level)) {
843 843
            // Calculation error
844 844
            _level->liftToTop(level);
845 845
          }
846 846
        } else {
847 847
          _level->deactivate(n);
848 848
        }
849 849

	
850 850
      }
851 851
    }
852 852

	
853 853
    /// \brief Runs the preflow algorithm.
854 854
    ///
855 855
    /// Runs the preflow algorithm.
856 856
    /// \note pf.run() is just a shortcut of the following code.
857 857
    /// \code
858 858
    ///   pf.init();
859 859
    ///   pf.startFirstPhase();
860 860
    ///   pf.startSecondPhase();
861 861
    /// \endcode
862 862
    void run() {
863 863
      init();
864 864
      startFirstPhase();
865 865
      startSecondPhase();
866 866
    }
867 867

	
868 868
    /// \brief Runs the preflow algorithm to compute the minimum cut.
869 869
    ///
870 870
    /// Runs the preflow algorithm to compute the minimum cut.
871 871
    /// \note pf.runMinCut() is just a shortcut of the following code.
872 872
    /// \code
873 873
    ///   pf.init();
874 874
    ///   pf.startFirstPhase();
875 875
    /// \endcode
876 876
    void runMinCut() {
877 877
      init();
878 878
      startFirstPhase();
879 879
    }
880 880

	
881 881
    /// @}
882 882

	
883 883
    /// \name Query Functions
884 884
    /// The results of the preflow algorithm can be obtained using these
885 885
    /// functions.\n
886 886
    /// Either one of the \ref run() "run*()" functions or one of the
887 887
    /// \ref startFirstPhase() "start*()" functions should be called
888 888
    /// before using them.
889 889

	
890 890
    ///@{
891 891

	
892 892
    /// \brief Returns the value of the maximum flow.
893 893
    ///
894 894
    /// Returns the value of the maximum flow by returning the excess
895 895
    /// of the target node. This value equals to the value of
896 896
    /// the maximum flow already after the first phase of the algorithm.
897 897
    ///
898 898
    /// \pre Either \ref run() or \ref init() must be called before
899 899
    /// using this function.
900
    Flow flowValue() const {
900
    Value flowValue() const {
901 901
      return (*_excess)[_target];
902 902
    }
903 903

	
904
    /// \brief Returns the flow on the given arc.
904
    /// \brief Returns the flow value on the given arc.
905 905
    ///
906
    /// Returns the flow on the given arc. This method can
906
    /// Returns the flow value on the given arc. This method can
907 907
    /// be called after the second phase of the algorithm.
908 908
    ///
909 909
    /// \pre Either \ref run() or \ref init() must be called before
910 910
    /// using this function.
911
    Flow flow(const Arc& arc) const {
911
    Value flow(const Arc& arc) const {
912 912
      return (*_flow)[arc];
913 913
    }
914 914

	
915 915
    /// \brief Returns a const reference to the flow map.
916 916
    ///
917 917
    /// Returns a const reference to the arc map storing the found flow.
918 918
    /// This method can be called after the second phase of the algorithm.
919 919
    ///
920 920
    /// \pre Either \ref run() or \ref init() must be called before
921 921
    /// using this function.
922 922
    const FlowMap& flowMap() const {
923 923
      return *_flow;
924 924
    }
925 925

	
926 926
    /// \brief Returns \c true when the node is on the source side of the
927 927
    /// minimum cut.
928 928
    ///
929 929
    /// Returns true when the node is on the source side of the found
930 930
    /// minimum cut. This method can be called both after running \ref
931 931
    /// startFirstPhase() and \ref startSecondPhase().
932 932
    ///
933 933
    /// \pre Either \ref run() or \ref init() must be called before
934 934
    /// using this function.
935 935
    bool minCut(const Node& node) const {
936 936
      return ((*_level)[node] == _level->maxLevel()) == _phase;
937 937
    }
938 938

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

	
961 961
    /// @}
962 962
  };
963 963
}
964 964

	
965 965
#endif
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
#include <sstream>
20 20
#include <lemon/lp_skeleton.h>
21 21
#include "test_tools.h"
22 22
#include <lemon/tolerance.h>
23 23

	
24
#ifdef HAVE_CONFIG_H
25 24
#include <lemon/config.h>
26
#endif
27 25

	
28 26
#ifdef LEMON_HAVE_GLPK
29 27
#include <lemon/glpk.h>
30 28
#endif
31 29

	
32 30
#ifdef LEMON_HAVE_CPLEX
33 31
#include <lemon/cplex.h>
34 32
#endif
35 33

	
36 34
#ifdef LEMON_HAVE_SOPLEX
37 35
#include <lemon/soplex.h>
38 36
#endif
39 37

	
40 38
#ifdef LEMON_HAVE_CLP
41 39
#include <lemon/clp.h>
42 40
#endif
43 41

	
44 42
using namespace lemon;
45 43

	
46 44
void lpTest(LpSolver& lp)
47 45
{
48 46

	
49 47
  typedef LpSolver LP;
50 48

	
51 49
  std::vector<LP::Col> x(10);
52 50
  //  for(int i=0;i<10;i++) x.push_back(lp.addCol());
53 51
  lp.addColSet(x);
54 52
  lp.colLowerBound(x,1);
55 53
  lp.colUpperBound(x,1);
56 54
  lp.colBounds(x,1,2);
57 55

	
58 56
  std::vector<LP::Col> y(10);
59 57
  lp.addColSet(y);
60 58

	
61 59
  lp.colLowerBound(y,1);
62 60
  lp.colUpperBound(y,1);
63 61
  lp.colBounds(y,1,2);
64 62

	
65 63
  std::map<int,LP::Col> z;
66 64

	
67 65
  z.insert(std::make_pair(12,INVALID));
68 66
  z.insert(std::make_pair(2,INVALID));
69 67
  z.insert(std::make_pair(7,INVALID));
70 68
  z.insert(std::make_pair(5,INVALID));
71 69

	
72 70
  lp.addColSet(z);
73 71

	
74 72
  lp.colLowerBound(z,1);
75 73
  lp.colUpperBound(z,1);
76 74
  lp.colBounds(z,1,2);
77 75

	
78 76
  {
79 77
    LP::Expr e,f,g;
80 78
    LP::Col p1,p2,p3,p4,p5;
81 79
    LP::Constr c;
82 80

	
83 81
    p1=lp.addCol();
84 82
    p2=lp.addCol();
85 83
    p3=lp.addCol();
86 84
    p4=lp.addCol();
87 85
    p5=lp.addCol();
88 86

	
89 87
    e[p1]=2;
90 88
    *e=12;
91 89
    e[p1]+=2;
92 90
    *e+=12;
93 91
    e[p1]-=2;
94 92
    *e-=12;
95 93

	
96 94
    e=2;
97 95
    e=2.2;
98 96
    e=p1;
99 97
    e=f;
100 98

	
101 99
    e+=2;
102 100
    e+=2.2;
103 101
    e+=p1;
104 102
    e+=f;
105 103

	
106 104
    e-=2;
107 105
    e-=2.2;
108 106
    e-=p1;
109 107
    e-=f;
110 108

	
111 109
    e*=2;
112 110
    e*=2.2;
113 111
    e/=2;
114 112
    e/=2.2;
115 113

	
116 114
    e=((p1+p2)+(p1-p2)+(p1+12)+(12+p1)+(p1-12)+(12-p1)+
117 115
       (f+12)+(12+f)+(p1+f)+(f+p1)+(f+g)+
118 116
       (f-12)+(12-f)+(p1-f)+(f-p1)+(f-g)+
119 117
       2.2*f+f*2.2+f/2.2+
120 118
       2*f+f*2+f/2+
121 119
       2.2*p1+p1*2.2+p1/2.2+
122 120
       2*p1+p1*2+p1/2
123 121
       );
124 122

	
125 123

	
126 124
    c = (e  <= f  );
127 125
    c = (e  <= 2.2);
128 126
    c = (e  <= 2  );
129 127
    c = (e  <= p1 );
130 128
    c = (2.2<= f  );
131 129
    c = (2  <= f  );
132 130
    c = (p1 <= f  );
133 131
    c = (p1 <= p2 );
134 132
    c = (p1 <= 2.2);
135 133
    c = (p1 <= 2  );
136 134
    c = (2.2<= p2 );
137 135
    c = (2  <= p2 );
138 136

	
139 137
    c = (e  >= f  );
140 138
    c = (e  >= 2.2);
141 139
    c = (e  >= 2  );
142 140
    c = (e  >= p1 );
143 141
    c = (2.2>= f  );
144 142
    c = (2  >= f  );
145 143
    c = (p1 >= f  );
146 144
    c = (p1 >= p2 );
147 145
    c = (p1 >= 2.2);
148 146
    c = (p1 >= 2  );
149 147
    c = (2.2>= p2 );
150 148
    c = (2  >= p2 );
151 149

	
152 150
    c = (e  == f  );
153 151
    c = (e  == 2.2);
154 152
    c = (e  == 2  );
155 153
    c = (e  == p1 );
156 154
    c = (2.2== f  );
157 155
    c = (2  == f  );
158 156
    c = (p1 == f  );
159 157
    //c = (p1 == p2 );
160 158
    c = (p1 == 2.2);
161 159
    c = (p1 == 2  );
162 160
    c = (2.2== p2 );
163 161
    c = (2  == p2 );
164 162

	
165 163
    c = ((2 <= e) <= 3);
166 164
    c = ((2 <= p1) <= 3);
167 165

	
168 166
    c = ((2 >= e) >= 3);
169 167
    c = ((2 >= p1) >= 3);
170 168

	
171 169
    e[x[3]]=2;
172 170
    e[x[3]]=4;
173 171
    e[x[3]]=1;
174 172
    *e=12;
175 173

	
176 174
    lp.addRow(-LP::INF,e,23);
177 175
    lp.addRow(-LP::INF,3.0*(x[1]+x[2]/2)-x[3],23);
178 176
    lp.addRow(-LP::INF,3.0*(x[1]+x[2]*2-5*x[3]+12-x[4]/3)+2*x[4]-4,23);
179 177

	
180 178
    lp.addRow(x[1]+x[3]<=x[5]-3);
181 179
    lp.addRow((-7<=x[1]+x[3]-12)<=3);
182 180
    lp.addRow(x[1]<=x[5]);
183 181

	
184 182
    std::ostringstream buf;
185 183

	
186 184

	
187 185
    e=((p1+p2)+(p1-0.99*p2));
188 186
    //e.prettyPrint(std::cout);
189 187
    //(e<=2).prettyPrint(std::cout);
190 188
    double tolerance=0.001;
191 189
    e.simplify(tolerance);
192 190
    buf << "Coeff. of p2 should be 0.01";
193 191
    check(e[p2]>0, buf.str());
194 192

	
195 193
    tolerance=0.02;
196 194
    e.simplify(tolerance);
197 195
    buf << "Coeff. of p2 should be 0";
198 196
    check(const_cast<const LpSolver::Expr&>(e)[p2]==0, buf.str());
199 197

	
200 198
    //Test for clone/new
201 199
    LP* lpnew = lp.newSolver();
202 200
    LP* lpclone = lp.cloneSolver();
203 201
    delete lpnew;
204 202
    delete lpclone;
205 203

	
206 204
  }
207 205

	
208 206
  {
209 207
    LP::DualExpr e,f,g;
210 208
    LP::Row p1 = INVALID, p2 = INVALID, p3 = INVALID,
211 209
      p4 = INVALID, p5 = INVALID;
212 210

	
213 211
    e[p1]=2;
214 212
    e[p1]+=2;
215 213
    e[p1]-=2;
216 214

	
217 215
    e=p1;
218 216
    e=f;
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
#include <iostream>
20 20
#include <fstream>
21
#include <limits>
21 22

	
22 23
#include <lemon/list_graph.h>
23 24
#include <lemon/lgf_reader.h>
24 25

	
25 26
#include <lemon/network_simplex.h>
26 27

	
27 28
#include <lemon/concepts/digraph.h>
28 29
#include <lemon/concept_check.h>
29 30

	
30 31
#include "test_tools.h"
31 32

	
32 33
using namespace lemon;
33 34

	
34 35
char test_lgf[] =
35 36
  "@nodes\n"
36
  "label  sup1 sup2 sup3 sup4 sup5\n"
37
  "    1    20   27    0   20   30\n"
38
  "    2    -4    0    0   -8   -3\n"
39
  "    3     0    0    0    0    0\n"
40
  "    4     0    0    0    0    0\n"
41
  "    5     9    0    0    6   11\n"
42
  "    6    -6    0    0   -5   -6\n"
43
  "    7     0    0    0    0    0\n"
44
  "    8     0    0    0    0    3\n"
45
  "    9     3    0    0    0    0\n"
46
  "   10    -2    0    0   -7   -2\n"
47
  "   11     0    0    0  -10    0\n"
48
  "   12   -20  -27    0  -30  -20\n"
49
  "\n"
37
  "label  sup1 sup2 sup3 sup4 sup5 sup6\n"
38
  "    1    20   27    0   30   20   30\n"
39
  "    2    -4    0    0    0   -8   -3\n"
40
  "    3     0    0    0    0    0    0\n"
41
  "    4     0    0    0    0    0    0\n"
42
  "    5     9    0    0    0    6   11\n"
43
  "    6    -6    0    0    0   -5   -6\n"
44
  "    7     0    0    0    0    0    0\n"
45
  "    8     0    0    0    0    0    3\n"
46
  "    9     3    0    0    0    0    0\n"
47
  "   10    -2    0    0    0   -7   -2\n"
48
  "   11     0    0    0    0  -10    0\n"
49
  "   12   -20  -27    0  -30  -30  -20\n"
50
  "\n"                
50 51
  "@arcs\n"
51
  "       cost  cap low1 low2\n"
52
  " 1  2    70   11    0    8\n"
53
  " 1  3   150    3    0    1\n"
54
  " 1  4    80   15    0    2\n"
55
  " 2  8    80   12    0    0\n"
56
  " 3  5   140    5    0    3\n"
57
  " 4  6    60   10    0    1\n"
58
  " 4  7    80    2    0    0\n"
59
  " 4  8   110    3    0    0\n"
60
  " 5  7    60   14    0    0\n"
61
  " 5 11   120   12    0    0\n"
62
  " 6  3     0    3    0    0\n"
63
  " 6  9   140    4    0    0\n"
64
  " 6 10    90    8    0    0\n"
65
  " 7  1    30    5    0    0\n"
66
  " 8 12    60   16    0    4\n"
67
  " 9 12    50    6    0    0\n"
68
  "10 12    70   13    0    5\n"
69
  "10  2   100    7    0    0\n"
70
  "10  7    60   10    0    0\n"
71
  "11 10    20   14    0    6\n"
72
  "12 11    30   10    0    0\n"
52
  "       cost  cap low1 low2 low3\n"
53
  " 1  2    70   11    0    8    8\n"
54
  " 1  3   150    3    0    1    0\n"
55
  " 1  4    80   15    0    2    2\n"
56
  " 2  8    80   12    0    0    0\n"
57
  " 3  5   140    5    0    3    1\n"
58
  " 4  6    60   10    0    1    0\n"
59
  " 4  7    80    2    0    0    0\n"
60
  " 4  8   110    3    0    0    0\n"
61
  " 5  7    60   14    0    0    0\n"
62
  " 5 11   120   12    0    0    0\n"
63
  " 6  3     0    3    0    0    0\n"
64
  " 6  9   140    4    0    0    0\n"
65
  " 6 10    90    8    0    0    0\n"
66
  " 7  1    30    5    0    0   -5\n"
67
  " 8 12    60   16    0    4    3\n"
68
  " 9 12    50    6    0    0    0\n"
69
  "10 12    70   13    0    5    2\n"
70
  "10  2   100    7    0    0    0\n"
71
  "10  7    60   10    0    0   -3\n"
72
  "11 10    20   14    0    6  -20\n"
73
  "12 11    30   10    0    0  -10\n"
73 74
  "\n"
74 75
  "@attributes\n"
75 76
  "source 1\n"
76 77
  "target 12\n";
77 78

	
78 79

	
79
enum ProblemType {
80
enum SupplyType {
80 81
  EQ,
81 82
  GEQ,
82 83
  LEQ
83 84
};
84 85

	
85 86
// Check the interface of an MCF algorithm
86
template <typename GR, typename Flow, typename Cost>
87
template <typename GR, typename Value, typename Cost>
87 88
class McfClassConcept
88 89
{
89 90
public:
90 91

	
91 92
  template <typename MCF>
92 93
  struct Constraints {
93 94
    void constraints() {
94 95
      checkConcept<concepts::Digraph, GR>();
95 96

	
96 97
      MCF mcf(g);
98
      const MCF& const_mcf = mcf;
97 99

	
98 100
      b = mcf.reset()
99 101
             .lowerMap(lower)
100 102
             .upperMap(upper)
101
             .capacityMap(upper)
102
             .boundMaps(lower, upper)
103 103
             .costMap(cost)
104 104
             .supplyMap(sup)
105 105
             .stSupply(n, n, k)
106
             .flowMap(flow)
107
             .potentialMap(pot)
108 106
             .run();
109
      
110
      const MCF& const_mcf = mcf;
111 107

	
112
      const typename MCF::FlowMap &fm = const_mcf.flowMap();
113
      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
114

	
115
      v = const_mcf.totalCost();
116
      double x = const_mcf.template totalCost<double>();
108
      c = const_mcf.totalCost();
109
      x = const_mcf.template totalCost<double>();
117 110
      v = const_mcf.flow(a);
118
      v = const_mcf.potential(n);
119

	
120
      ignore_unused_variable_warning(fm);
121
      ignore_unused_variable_warning(pm);
122
      ignore_unused_variable_warning(x);
111
      c = const_mcf.potential(n);
112
      const_mcf.flowMap(fm);
113
      const_mcf.potentialMap(pm);
123 114
    }
124 115

	
125 116
    typedef typename GR::Node Node;
126 117
    typedef typename GR::Arc Arc;
127
    typedef concepts::ReadMap<Node, Flow> NM;
128
    typedef concepts::ReadMap<Arc, Flow> FAM;
118
    typedef concepts::ReadMap<Node, Value> NM;
119
    typedef concepts::ReadMap<Arc, Value> VAM;
129 120
    typedef concepts::ReadMap<Arc, Cost> CAM;
121
    typedef concepts::WriteMap<Arc, Value> FlowMap;
122
    typedef concepts::WriteMap<Node, Cost> PotMap;
130 123

	
131 124
    const GR &g;
132
    const FAM &lower;
133
    const FAM &upper;
125
    const VAM &lower;
126
    const VAM &upper;
134 127
    const CAM &cost;
135 128
    const NM &sup;
136 129
    const Node &n;
137 130
    const Arc &a;
138
    const Flow &k;
139
    Flow v;
131
    const Value &k;
132
    FlowMap fm;
133
    PotMap pm;
140 134
    bool b;
141

	
142
    typename MCF::FlowMap &flow;
143
    typename MCF::PotentialMap &pot;
135
    double x;
136
    typename MCF::Value v;
137
    typename MCF::Cost c;
144 138
  };
145 139

	
146 140
};
147 141

	
148 142

	
149 143
// Check the feasibility of the given flow (primal soluiton)
150 144
template < typename GR, typename LM, typename UM,
151 145
           typename SM, typename FM >
152 146
bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
153 147
                const SM& supply, const FM& flow,
154
                ProblemType type = EQ )
148
                SupplyType type = EQ )
155 149
{
156 150
  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
157 151

	
158 152
  for (ArcIt e(gr); e != INVALID; ++e) {
159 153
    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
160 154
  }
161 155

	
162 156
  for (NodeIt n(gr); n != INVALID; ++n) {
163 157
    typename SM::Value sum = 0;
164 158
    for (OutArcIt e(gr, n); e != INVALID; ++e)
165 159
      sum += flow[e];
166 160
    for (InArcIt e(gr, n); e != INVALID; ++e)
167 161
      sum -= flow[e];
168 162
    bool b = (type ==  EQ && sum == supply[n]) ||
169 163
             (type == GEQ && sum >= supply[n]) ||
170 164
             (type == LEQ && sum <= supply[n]);
171 165
    if (!b) return false;
172 166
  }
173 167

	
174 168
  return true;
175 169
}
176 170

	
177 171
// Check the feasibility of the given potentials (dual soluiton)
178 172
// using the "Complementary Slackness" optimality condition
179 173
template < typename GR, typename LM, typename UM,
180 174
           typename CM, typename SM, typename FM, typename PM >
181 175
bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
182 176
                     const CM& cost, const SM& supply, const FM& flow, 
183 177
                     const PM& pi )
184 178
{
185 179
  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
186 180

	
187 181
  bool opt = true;
188 182
  for (ArcIt e(gr); opt && e != INVALID; ++e) {
189 183
    typename CM::Value red_cost =
190 184
      cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
191 185
    opt = red_cost == 0 ||
192 186
          (red_cost > 0 && flow[e] == lower[e]) ||
193 187
          (red_cost < 0 && flow[e] == upper[e]);
194 188
  }
195 189
  
196 190
  for (NodeIt n(gr); opt && n != INVALID; ++n) {
197 191
    typename SM::Value sum = 0;
198 192
    for (OutArcIt e(gr, n); e != INVALID; ++e)
199 193
      sum += flow[e];
200 194
    for (InArcIt e(gr, n); e != INVALID; ++e)
201 195
      sum -= flow[e];
202 196
    opt = (sum == supply[n]) || (pi[n] == 0);
203 197
  }
204 198
  
205 199
  return opt;
206 200
}
207 201

	
208 202
// Run a minimum cost flow algorithm and check the results
209 203
template < typename MCF, typename GR,
210 204
           typename LM, typename UM,
211
           typename CM, typename SM >
212
void checkMcf( const MCF& mcf, bool mcf_result,
205
           typename CM, typename SM,
206
           typename PT >
207
void checkMcf( const MCF& mcf, PT mcf_result,
213 208
               const GR& gr, const LM& lower, const UM& upper,
214 209
               const CM& cost, const SM& supply,
215
               bool result, typename CM::Value total,
210
               PT result, bool optimal, typename CM::Value total,
216 211
               const std::string &test_id = "",
217
               ProblemType type = EQ )
212
               SupplyType type = EQ )
218 213
{
219 214
  check(mcf_result == result, "Wrong result " + test_id);
220
  if (result) {
221
    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
215
  if (optimal) {
216
    typename GR::template ArcMap<typename SM::Value> flow(gr);
217
    typename GR::template NodeMap<typename CM::Value> pi(gr);
218
    mcf.flowMap(flow);
219
    mcf.potentialMap(pi);
220
    check(checkFlow(gr, lower, upper, supply, flow, type),
222 221
          "The flow is not feasible " + test_id);
223 222
    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
224
    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
225
                         mcf.potentialMap()),
223
    check(checkPotential(gr, lower, upper, cost, supply, flow, pi),
226 224
          "Wrong potentials " + test_id);
227 225
  }
228 226
}
229 227

	
230 228
int main()
231 229
{
232 230
  // Check the interfaces
233 231
  {
234
    typedef int Flow;
235
    typedef int Cost;
236 232
    typedef concepts::Digraph GR;
237
    checkConcept< McfClassConcept<GR, Flow, Cost>,
238
                  NetworkSimplex<GR, Flow, Cost> >();
233
    checkConcept< McfClassConcept<GR, int, int>,
234
                  NetworkSimplex<GR> >();
235
    checkConcept< McfClassConcept<GR, double, double>,
236
                  NetworkSimplex<GR, double> >();
237
    checkConcept< McfClassConcept<GR, int, double>,
238
                  NetworkSimplex<GR, int, double> >();
239 239
  }
240 240

	
241 241
  // Run various MCF tests
242 242
  typedef ListDigraph Digraph;
243 243
  DIGRAPH_TYPEDEFS(ListDigraph);
244 244

	
245 245
  // Read the test digraph
246 246
  Digraph gr;
247
  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
248
  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
247
  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
248
  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
249 249
  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
250 250
  Node v, w;
251 251

	
252 252
  std::istringstream input(test_lgf);
253 253
  DigraphReader<Digraph>(gr, input)
254 254
    .arcMap("cost", c)
255 255
    .arcMap("cap", u)
256 256
    .arcMap("low1", l1)
257 257
    .arcMap("low2", l2)
258
    .arcMap("low3", l3)
258 259
    .nodeMap("sup1", s1)
259 260
    .nodeMap("sup2", s2)
260 261
    .nodeMap("sup3", s3)
261 262
    .nodeMap("sup4", s4)
262 263
    .nodeMap("sup5", s5)
264
    .nodeMap("sup6", s6)
263 265
    .node("source", v)
264 266
    .node("target", w)
265 267
    .run();
268
  
269
  // Build a test digraph for testing negative costs
270
  Digraph ngr;
271
  Node n1 = ngr.addNode();
272
  Node n2 = ngr.addNode();
273
  Node n3 = ngr.addNode();
274
  Node n4 = ngr.addNode();
275
  Node n5 = ngr.addNode();
276
  Node n6 = ngr.addNode();
277
  Node n7 = ngr.addNode();
278
  
279
  Arc a1 = ngr.addArc(n1, n2);
280
  Arc a2 = ngr.addArc(n1, n3);
281
  Arc a3 = ngr.addArc(n2, n4);
282
  Arc a4 = ngr.addArc(n3, n4);
283
  Arc a5 = ngr.addArc(n3, n2);
284
  Arc a6 = ngr.addArc(n5, n3);
285
  Arc a7 = ngr.addArc(n5, n6);
286
  Arc a8 = ngr.addArc(n6, n7);
287
  Arc a9 = ngr.addArc(n7, n5);
288
  
289
  Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
290
  ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
291
  Digraph::NodeMap<int> ns(ngr, 0);
292
  
293
  nl2[a7] =  1000;
294
  nl2[a8] = -1000;
295
  
296
  ns[n1] =  100;
297
  ns[n4] = -100;
298
  
299
  nc[a1] =  100;
300
  nc[a2] =   30;
301
  nc[a3] =   20;
302
  nc[a4] =   80;
303
  nc[a5] =   50;
304
  nc[a6] =   10;
305
  nc[a7] =   80;
306
  nc[a8] =   30;
307
  nc[a9] = -120;
266 308

	
267 309
  // A. Test NetworkSimplex with the default pivot rule
268 310
  {
269 311
    NetworkSimplex<Digraph> mcf(gr);
270 312

	
271 313
    // Check the equality form
272 314
    mcf.upperMap(u).costMap(c);
273 315
    checkMcf(mcf, mcf.supplyMap(s1).run(),
274
             gr, l1, u, c, s1, true,  5240, "#A1");
316
             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
275 317
    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
276
             gr, l1, u, c, s2, true,  7620, "#A2");
318
             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
277 319
    mcf.lowerMap(l2);
278 320
    checkMcf(mcf, mcf.supplyMap(s1).run(),
279
             gr, l2, u, c, s1, true,  5970, "#A3");
321
             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
280 322
    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
281
             gr, l2, u, c, s2, true,  8010, "#A4");
323
             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
282 324
    mcf.reset();
283 325
    checkMcf(mcf, mcf.supplyMap(s1).run(),
284
             gr, l1, cu, cc, s1, true,  74, "#A5");
326
             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
285 327
    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
286
             gr, l2, cu, cc, s2, true,  94, "#A6");
328
             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
287 329
    mcf.reset();
288 330
    checkMcf(mcf, mcf.run(),
289
             gr, l1, cu, cc, s3, true,   0, "#A7");
290
    checkMcf(mcf, mcf.boundMaps(l2, u).run(),
291
             gr, l2, u, cc, s3, false,   0, "#A8");
331
             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
332
    checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
333
             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
334
    mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
335
    checkMcf(mcf, mcf.run(),
336
             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
292 337

	
293 338
    // Check the GEQ form
294
    mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
339
    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
295 340
    checkMcf(mcf, mcf.run(),
296
             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
297
    mcf.problemType(mcf.GEQ);
341
             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
342
    mcf.supplyType(mcf.GEQ);
298 343
    checkMcf(mcf, mcf.lowerMap(l2).run(),
299
             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
300
    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
344
             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
345
    mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
301 346
    checkMcf(mcf, mcf.run(),
302
             gr, l2, u, c, s5, false,    0, "#A11", GEQ);
347
             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
303 348

	
304 349
    // Check the LEQ form
305
    mcf.reset().problemType(mcf.LEQ);
306
    mcf.upperMap(u).costMap(c).supplyMap(s5);
350
    mcf.reset().supplyType(mcf.LEQ);
351
    mcf.upperMap(u).costMap(c).supplyMap(s6);
307 352
    checkMcf(mcf, mcf.run(),
308
             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
353
             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
309 354
    checkMcf(mcf, mcf.lowerMap(l2).run(),
310
             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
311
    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
355
             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
356
    mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
312 357
    checkMcf(mcf, mcf.run(),
313
             gr, l2, u, c, s4, false,    0, "#A14", LEQ);
358
             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
359

	
360
    // Check negative costs
361
    NetworkSimplex<Digraph> nmcf(ngr);
362
    nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
363
    checkMcf(nmcf, nmcf.run(),
364
      ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
365
    checkMcf(nmcf, nmcf.upperMap(nu2).run(),
366
      ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
367
    nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
368
    checkMcf(nmcf, nmcf.run(),
369
      ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
314 370
  }
315 371

	
316 372
  // B. Test NetworkSimplex with each pivot rule
317 373
  {
318 374
    NetworkSimplex<Digraph> mcf(gr);
319
    mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
375
    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
320 376

	
321 377
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
322
             gr, l2, u, c, s1, true,  5970, "#B1");
378
             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
323 379
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
324
             gr, l2, u, c, s1, true,  5970, "#B2");
380
             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
325 381
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
326
             gr, l2, u, c, s1, true,  5970, "#B3");
382
             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
327 383
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
328
             gr, l2, u, c, s1, true,  5970, "#B4");
384
             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
329 385
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
330
             gr, l2, u, c, s1, true,  5970, "#B5");
386
             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
331 387
  }
332 388

	
333 389
  return 0;
334 390
}
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
#include "test_tools.h"
20 20

	
21
#ifdef HAVE_CONFIG_H
22 21
#include <lemon/config.h>
23
#endif
24 22

	
25 23
#ifdef LEMON_HAVE_CPLEX
26 24
#include <lemon/cplex.h>
27 25
#endif
28 26

	
29 27
#ifdef LEMON_HAVE_GLPK
30 28
#include <lemon/glpk.h>
31 29
#endif
32 30

	
33 31
#ifdef LEMON_HAVE_CBC
34 32
#include <lemon/cbc.h>
35 33
#endif
36 34

	
37 35

	
38 36
using namespace lemon;
39 37

	
40 38
void solveAndCheck(MipSolver& mip, MipSolver::ProblemType stat,
41 39
                   double exp_opt) {
42 40
  using std::string;
43 41

	
44 42
  mip.solve();
45 43
  //int decimal,sign;
46 44
  std::ostringstream buf;
47 45
  buf << "Type should be: " << int(stat)<<" and it is "<<int(mip.type());
48 46

	
49 47

	
50 48
  //  itoa(stat,buf1, 10);
51 49
  check(mip.type()==stat, buf.str());
52 50

	
53 51
  if (stat ==  MipSolver::OPTIMAL) {
54 52
    std::ostringstream sbuf;
55 53
    buf << "Wrong optimal value: the right optimum is " << exp_opt;
56 54
    check(std::abs(mip.solValue()-exp_opt) < 1e-3, sbuf.str());
57 55
    //+ecvt(exp_opt,2)
58 56
  }
59 57
}
60 58

	
61 59
void aTest(MipSolver& mip)
62 60
{
63 61
  //The following example is very simple
64 62

	
65 63

	
66 64
  typedef MipSolver::Row Row;
67 65
  typedef MipSolver::Col Col;
68 66

	
69 67

	
70 68
  Col x1 = mip.addCol();
71 69
  Col x2 = mip.addCol();
72 70

	
73 71

	
74 72
  //Objective function
75 73
  mip.obj(x1);
76 74

	
77 75
  mip.max();
78 76

	
79 77
  //Unconstrained optimization
80 78
  mip.solve();
81 79
  //Check it out!
82 80

	
83 81
  //Constraints
84 82
  mip.addRow(2 * x1 + x2 <= 2);
85 83
  Row y2 = mip.addRow(x1 - 2 * x2 <= 0);
86 84

	
87 85
  //Nonnegativity of the variable x1
88 86
  mip.colLowerBound(x1, 0);
89 87

	
90 88

	
91 89
  //Maximization of x1
92 90
  //over the triangle with vertices (0,0),(4/5,2/5),(0,2)
93 91
  double expected_opt=4.0/5.0;
94 92
  solveAndCheck(mip, MipSolver::OPTIMAL, expected_opt);
95 93

	
96 94

	
97 95
  //Restrict x2 to integer
98 96
  mip.colType(x2,MipSolver::INTEGER);
99 97
  expected_opt=1.0/2.0;
100 98
  solveAndCheck(mip, MipSolver::OPTIMAL, expected_opt);
101 99

	
102 100

	
103 101
  //Restrict both to integer
104 102
  mip.colType(x1,MipSolver::INTEGER);
105 103
  expected_opt=0;
106 104
  solveAndCheck(mip, MipSolver::OPTIMAL, expected_opt);
107 105

	
108 106
  //Erase a variable
109 107
  mip.erase(x2);
110 108
  mip.rowUpperBound(y2, 8);
111 109
  expected_opt=1;
112 110
  solveAndCheck(mip, MipSolver::OPTIMAL, expected_opt);
113 111

	
114 112
}
115 113

	
116 114

	
117 115
template<class MIP>
118 116
void cloneTest()
119 117
{
120 118

	
121 119
  MIP* mip = new MIP();
122 120
  MIP* mipnew = mip->newSolver();
123 121
  MIP* mipclone = mip->cloneSolver();
124 122
  delete mip;
125 123
  delete mipnew;
126 124
  delete mipclone;
127 125
}
128 126

	
129 127
int main()
130 128
{
131 129

	
132 130
#ifdef LEMON_HAVE_GLPK
133 131
  {
134 132
    GlpkMip mip1;
135 133
    aTest(mip1);
136 134
    cloneTest<GlpkMip>();
137 135
  }
138 136
#endif
139 137

	
140 138
#ifdef LEMON_HAVE_CPLEX
141 139
  try {
142 140
    CplexMip mip2;
143 141
    aTest(mip2);
144 142
    cloneTest<CplexMip>();
145 143
  } catch (CplexEnv::LicenseError& error) {
146 144
    check(false, error.what());
147 145
  }
148 146
#endif
149 147

	
150 148
#ifdef LEMON_HAVE_CBC
151 149
  {
152 150
    CbcMip mip1;
153 151
    aTest(mip1);
154 152
    cloneTest<CbcMip>();
155 153
  }
156 154
#endif
157 155

	
158 156
  return 0;
159 157

	
160 158
}
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
///\ingroup tools
20 20
///\file
21 21
///\brief DIMACS problem solver.
22 22
///
23 23
/// This program solves various problems given in DIMACS format.
24 24
///
25 25
/// See
26 26
/// \code
27 27
///   dimacs-solver --help
28 28
/// \endcode
29 29
/// for more info on usage.
30 30

	
31 31
#include <iostream>
32 32
#include <fstream>
33 33
#include <cstring>
34 34

	
35 35
#include <lemon/smart_graph.h>
36 36
#include <lemon/dimacs.h>
37 37
#include <lemon/lgf_writer.h>
38 38
#include <lemon/time_measure.h>
39 39

	
40 40
#include <lemon/arg_parser.h>
41 41
#include <lemon/error.h>
42 42

	
43 43
#include <lemon/dijkstra.h>
44 44
#include <lemon/preflow.h>
45 45
#include <lemon/matching.h>
46 46
#include <lemon/network_simplex.h>
47 47

	
48 48
using namespace lemon;
49 49
typedef SmartDigraph Digraph;
50 50
DIGRAPH_TYPEDEFS(Digraph);
51 51
typedef SmartGraph Graph;
52 52

	
53 53
template<class Value>
54 54
void solve_sp(ArgParser &ap, std::istream &is, std::ostream &,
55 55
              DimacsDescriptor &desc)
56 56
{
57 57
  bool report = !ap.given("q");
58 58
  Digraph g;
59 59
  Node s;
60 60
  Digraph::ArcMap<Value> len(g);
61 61
  Timer t;
62 62
  t.restart();
63 63
  readDimacsSp(is, g, len, s, desc);
64 64
  if(report) std::cerr << "Read the file: " << t << '\n';
65 65
  t.restart();
66 66
  Dijkstra<Digraph, Digraph::ArcMap<Value> > dij(g,len);
67 67
  if(report) std::cerr << "Setup Dijkstra class: " << t << '\n';
68 68
  t.restart();
69 69
  dij.run(s);
70 70
  if(report) std::cerr << "Run Dijkstra: " << t << '\n';
71 71
}
72 72

	
73 73
template<class Value>
74 74
void solve_max(ArgParser &ap, std::istream &is, std::ostream &,
75 75
               Value infty, DimacsDescriptor &desc)
76 76
{
77 77
  bool report = !ap.given("q");
78 78
  Digraph g;
79 79
  Node s,t;
80 80
  Digraph::ArcMap<Value> cap(g);
81 81
  Timer ti;
82 82
  ti.restart();
83 83
  readDimacsMax(is, g, cap, s, t, infty, desc);
84 84
  if(report) std::cerr << "Read the file: " << ti << '\n';
85 85
  ti.restart();
86 86
  Preflow<Digraph, Digraph::ArcMap<Value> > pre(g,cap,s,t);
87 87
  if(report) std::cerr << "Setup Preflow class: " << ti << '\n';
88 88
  ti.restart();
89 89
  pre.run();
90 90
  if(report) std::cerr << "Run Preflow: " << ti << '\n';
91 91
  if(report) std::cerr << "\nMax flow value: " << pre.flowValue() << '\n';  
92 92
}
93 93

	
94 94
template<class Value>
95 95
void solve_min(ArgParser &ap, std::istream &is, std::ostream &,
96 96
               Value infty, DimacsDescriptor &desc)
97 97
{
98 98
  bool report = !ap.given("q");
99 99
  Digraph g;
100 100
  Digraph::ArcMap<Value> lower(g), cap(g), cost(g);
101 101
  Digraph::NodeMap<Value> sup(g);
102 102
  Timer ti;
103 103

	
104 104
  ti.restart();
105 105
  readDimacsMin(is, g, lower, cap, cost, sup, infty, desc);
106 106
  ti.stop();
107 107
  Value sum_sup = 0;
108 108
  for (Digraph::NodeIt n(g); n != INVALID; ++n) {
109 109
    sum_sup += sup[n];
110 110
  }
111 111
  if (report) {
112 112
    std::cerr << "Sum of supply values: " << sum_sup << "\n";
113 113
    if (sum_sup <= 0)
114 114
      std::cerr << "GEQ supply contraints are used for NetworkSimplex\n\n";
115 115
    else
116 116
      std::cerr << "LEQ supply contraints are used for NetworkSimplex\n\n";
117 117
  }
118 118
  if (report) std::cerr << "Read the file: " << ti << '\n';
119 119

	
120 120
  ti.restart();
121 121
  NetworkSimplex<Digraph, Value> ns(g);
122
  ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
123
  if (sum_sup > 0) ns.problemType(ns.LEQ);
122
  ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup);
123
  if (sum_sup > 0) ns.supplyType(ns.LEQ);
124 124
  if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';
125 125
  ti.restart();
126 126
  bool res = ns.run();
127 127
  if (report) {
128 128
    std::cerr << "Run NetworkSimplex: " << ti << "\n\n";
129 129
    std::cerr << "Feasible flow: " << (res ? "found" : "not found") << '\n';
130 130
    if (res) std::cerr << "Min flow cost: " << ns.totalCost() << '\n';
131 131
  }
132 132
}
133 133

	
134 134
void solve_mat(ArgParser &ap, std::istream &is, std::ostream &,
135 135
              DimacsDescriptor &desc)
136 136
{
137 137
  bool report = !ap.given("q");
138 138
  Graph g;
139 139
  Timer ti;
140 140
  ti.restart();
141 141
  readDimacsMat(is, g, desc);
142 142
  if(report) std::cerr << "Read the file: " << ti << '\n';
143 143
  ti.restart();
144 144
  MaxMatching<Graph> mat(g);
145 145
  if(report) std::cerr << "Setup MaxMatching class: " << ti << '\n';
146 146
  ti.restart();
147 147
  mat.run();
148 148
  if(report) std::cerr << "Run MaxMatching: " << ti << '\n';
149 149
  if(report) std::cerr << "\nCardinality of max matching: "
150 150
                       << mat.matchingSize() << '\n';  
151 151
}
152 152

	
153 153

	
154 154
template<class Value>
155 155
void solve(ArgParser &ap, std::istream &is, std::ostream &os,
156 156
           DimacsDescriptor &desc)
157 157
{
158 158
  std::stringstream iss(static_cast<std::string>(ap["infcap"]));
159 159
  Value infty;
160 160
  iss >> infty;
161 161
  if(iss.fail())
162 162
    {
163 163
      std::cerr << "Cannot interpret '"
164 164
                << static_cast<std::string>(ap["infcap"]) << "' as infinite"
165 165
                << std::endl;
166 166
      exit(1);
167 167
    }
168 168
  
169 169
  switch(desc.type)
170 170
    {
171 171
    case DimacsDescriptor::MIN:
172 172
      solve_min<Value>(ap,is,os,infty,desc);
173 173
      break;
174 174
    case DimacsDescriptor::MAX:
175 175
      solve_max<Value>(ap,is,os,infty,desc);
176 176
      break;
177 177
    case DimacsDescriptor::SP:
178 178
      solve_sp<Value>(ap,is,os,desc);
179 179
      break;
180 180
    case DimacsDescriptor::MAT:
181 181
      solve_mat(ap,is,os,desc);
182 182
      break;
183 183
    default:
184 184
      break;
185 185
    }
186 186
}
187 187

	
188 188
int main(int argc, const char *argv[]) {
189 189
  typedef SmartDigraph Digraph;
190 190

	
191 191
  typedef Digraph::Arc Arc;
192 192

	
193 193
  std::string inputName;
194 194
  std::string outputName;
195 195

	
196 196
  ArgParser ap(argc, argv);
197 197
  ap.other("[INFILE [OUTFILE]]",
198 198
           "If either the INFILE or OUTFILE file is missing the standard\n"
199 199
           "     input/output will be used instead.")
200 200
    .boolOption("q", "Do not print any report")
201 201
    .boolOption("int","Use 'int' for capacities, costs etc. (default)")
202 202
    .optionGroup("datatype","int")
203 203
#ifdef LEMON_HAVE_LONG_LONG
204 204
    .boolOption("long","Use 'long long' for capacities, costs etc.")
205 205
    .optionGroup("datatype","long")
206 206
#endif
207 207
    .boolOption("double","Use 'double' for capacities, costs etc.")
208 208
    .optionGroup("datatype","double")
209 209
    .boolOption("ldouble","Use 'long double' for capacities, costs etc.")
210 210
    .optionGroup("datatype","ldouble")
211 211
    .onlyOneGroup("datatype")
212 212
    .stringOption("infcap","Value used for 'very high' capacities","0")
213 213
    .run();
214 214

	
215 215
  std::ifstream input;
216 216
  std::ofstream output;
217 217

	
218 218
  switch(ap.files().size())
219 219
    {
220 220
    case 2:
221 221
      output.open(ap.files()[1].c_str());
222 222
      if (!output) {
223 223
        throw IoError("Cannot open the file for writing", ap.files()[1]);
224 224
      }
225 225
    case 1:
226 226
      input.open(ap.files()[0].c_str());
227 227
      if (!input) {
228 228
        throw IoError("File cannot be found", ap.files()[0]);
229 229
      }
230 230
    case 0:
231 231
      break;
232 232
    default:
233 233
      std::cerr << ap.commandName() << ": too many arguments\n";
234 234
      return 1;
235 235
    }
236 236
  std::istream& is = (ap.files().size()<1 ? std::cin : input);
237 237
  std::ostream& os = (ap.files().size()<2 ? std::cout : output);
238 238

	
239 239
  DimacsDescriptor desc = dimacsType(is);
240 240
  
241 241
  if(!ap.given("q"))
242 242
    {
243 243
      std::cout << "Problem type: ";
244 244
      switch(desc.type)
245 245
        {
246 246
        case DimacsDescriptor::MIN:
247 247
          std::cout << "min";
248 248
          break;
249 249
        case DimacsDescriptor::MAX:
250 250
          std::cout << "max";
251 251
          break;
252 252
        case DimacsDescriptor::SP:
253 253
          std::cout << "sp";
254 254
        case DimacsDescriptor::MAT:
255 255
          std::cout << "mat";
256 256
          break;
257 257
        default:
258 258
          exit(1);
259 259
          break;
260 260
        }
261 261
      std::cout << "\nNum of nodes: " << desc.nodeNum;
262 262
      std::cout << "\nNum of arcs:  " << desc.edgeNum;
263 263
      std::cout << "\n\n";
264 264
    }
265 265
    
266 266
  if(ap.given("double"))
267 267
    solve<double>(ap,is,os,desc);
268 268
  else if(ap.given("ldouble"))
269 269
    solve<long double>(ap,is,os,desc);
270 270
#ifdef LEMON_HAVE_LONG_LONG
271 271
  else if(ap.given("long"))
272 272
    solve<long long>(ap,is,os,desc);
273 273
#endif
274 274
  else solve<int>(ap,is,os,desc);
275 275

	
276 276
  return 0;
277 277
}
0 comments (0 inline)