↑ 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
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
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(
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
... ...
@@ -331,118 +331,118 @@
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
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 \
Ignore white space 6 line context
... ...
@@ -43,89 +43,89 @@
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$
... ...
@@ -166,83 +166,83 @@
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
... ...
@@ -509,110 +509,110 @@
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);
... ...
@@ -640,55 +640,55 @@
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

	
... ...
@@ -729,66 +729,66 @@
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

	
Ignore white space 6 line context
... ...
@@ -9,276 +9,280 @@
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) :
... ...
@@ -638,764 +642,574 @@
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];
... ...
@@ -1505,82 +1319,96 @@
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 6 line context
... ...
@@ -25,154 +25,154 @@
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() {
... ...
@@ -449,243 +449,243 @@
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;
... ...
@@ -757,77 +757,77 @@
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;
... ...
@@ -876,60 +876,60 @@
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 {
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

	
Ignore white space 48 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
... ...
@@ -187,148 +181,210 @@
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());
Ignore white space 6 line context
... ...
@@ -98,50 +98,50 @@
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();
0 comments (0 inline)