0
15
0
3
2
42
22
10
10
19
19
1
1
344
516
21
21
160
104
... | ... |
@@ -8,37 +8,33 @@ |
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}) |
... | ... |
@@ -2,29 +2,30 @@ |
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 = |
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 |
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) |
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 |
|
|
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) |
... | ... |
@@ -343,35 +343,35 @@ |
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 |
|
|
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} |
|
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, |
... | ... |
@@ -395,42 +395,42 @@ |
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} |
|
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 |
|
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 |
... | ... |
@@ -6,24 +6,25 @@ |
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 |
nodist_lemon_HEADERS = lemon/config.h |
|
18 | 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) \ |
... | ... |
@@ -48,29 +49,29 @@ |
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 \ |
... | ... |
@@ -55,33 +55,33 @@ |
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< |
|
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 |
/// |
... | ... |
@@ -95,25 +95,25 @@ |
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< |
|
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 |
... | ... |
@@ -178,26 +178,26 @@ |
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; |
... | ... |
@@ -212,25 +212,25 @@ |
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< |
|
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 |
///@{ |
... | ... |
@@ -521,25 +521,25 @@ |
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 |
|
|
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])) |
... | ... |
@@ -554,53 +554,53 @@ |
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 |
|
|
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 |
|
|
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 |
|
|
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 { |
... | ... |
@@ -652,31 +652,31 @@ |
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 |
|
|
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 |
} |
... | ... |
@@ -741,43 +741,43 @@ |
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 |
|
|
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]; |
... | ... |
@@ -13,25 +13,25 @@ |
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/ |
|
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 { |
... | ... |
@@ -21,87 +21,147 @@ |
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 |
|
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 |
|
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 |
|
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 |
|
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 |
|
82 |
/// \brief Problem type constants for the \c run() function. |
|
92 | 83 |
/// |
93 |
/// Enum type |
|
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 |
... | ... |
@@ -122,151 +182,95 @@ |
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< |
|
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 |
|
|
226 |
IntArcMap _arc_id; |
|
234 | 227 |
IntVector _source; |
235 | 228 |
IntVector _target; |
236 | 229 |
|
237 | 230 |
// Node and arc data |
238 |
|
|
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 |
|
|
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; |
... | ... |
@@ -650,667 +654,475 @@ |
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"); |
|
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; |
|
682 |
|
|
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; |
|
675 | 715 |
} |
676 | 716 |
|
677 |
/// Destructor. |
|
678 |
~NetworkSimplex() { |
|
679 |
if (_local_flow) delete _flow_map; |
|
680 |
if (_local_potential) delete _potential_map; |
|
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 |
/// |
|
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 |
|
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 |
|
|
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 |
|
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 |
|
|
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 |
|
|
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 |
|
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 |
|
|
807 |
template<typename SupplyMap> |
|
808 |
NetworkSimplex& supplyMap(const SupplyMap& map) { |
|
809 | 809 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
810 |
|
|
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 |
|
|
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 |
|
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 |
|
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; |
|
850 |
return *this; |
|
851 |
} |
|
852 |
|
|
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 = ↦ |
|
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 |
|
|
851 |
NetworkSimplex& supplyType(SupplyType supply_type) { |
|
852 |
_stype = supply_type; |
|
885 | 853 |
return *this; |
886 | 854 |
} |
887 | 855 |
|
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 |
|
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. |
|
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 |
|
|
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 |
|
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). |
|
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. |
|
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 |
/// |
|
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 |
|
|
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 |
|
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 |
|
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 |
/// |
|
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 |
|
|
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 |
|
|
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]; |
|
1109 | 1077 |
} |
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]; |
|
1116 |
} |
|
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] = |
|
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 && |
|
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] = - |
|
1112 |
_pi[u] = -ART_COST + _pi[_root]; |
|
1301 | 1113 |
} else { |
1302 | 1114 |
_flow[e] = -_supply[u]; |
1303 | 1115 |
_forward[u] = false; |
1304 |
_pi[u] = |
|
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]) { |
... | ... |
@@ -1327,63 +1139,65 @@ |
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 |
|
|
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] ? |
|
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] ? |
|
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 |
|
|
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]] = |
... | ... |
@@ -1517,70 +1331,84 @@ |
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 |
|
|
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 |
|
1357 |
return INFEASIBLE; // avoid warning |
|
1544 | 1358 |
} |
1545 | 1359 |
|
1546 | 1360 |
template <typename PivotRuleImpl> |
1547 |
|
|
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 | 1375 |
|
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 |
|
|
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; |
|
1566 | 1380 |
} |
1567 |
} else { |
|
1568 |
for (int i = 0; i != _arc_num; ++i) { |
|
1569 |
_flow_map->set(_arc_ref[i], _flow[i]); |
|
1570 | 1381 |
} |
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; |
|
1571 | 1385 |
} |
1572 |
// Copy potential values to _potential_map |
|
1573 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
|
1574 |
|
|
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 |
|
|
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 |
... | ... |
@@ -37,31 +37,31 @@ |
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 |
|
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< |
|
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 |
/// |
... | ... |
@@ -75,25 +75,25 @@ |
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< |
|
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 |
... | ... |
@@ -116,25 +116,25 @@ |
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:: |
|
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 |
|
... | ... |
@@ -142,25 +142,25 @@ |
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< |
|
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); |
... | ... |
@@ -461,25 +461,25 @@ |
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 |
|
|
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 |
|
... | ... |
@@ -510,37 +510,37 @@ |
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 |
|
|
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 |
|
|
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 |
} |
... | ... |
@@ -555,52 +555,52 @@ |
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 |
|
|
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 |
|
|
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 |
|
|
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 { |
... | ... |
@@ -628,52 +628,52 @@ |
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 |
|
|
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 |
|
|
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 |
|
|
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 { |
... | ... |
@@ -769,53 +769,53 @@ |
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 |
|
|
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 |
|
|
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 |
|
|
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 { |
... | ... |
@@ -888,36 +888,36 @@ |
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 |
|
|
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 |
|
|
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; |
... | ... |
@@ -12,27 +12,25 @@ |
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 |
... | ... |
@@ -9,158 +9,152 @@ |
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 |
" |
|
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" |
|
49 | 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 |
|
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 |
|
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 | 107 |
|
110 |
const MCF& const_mcf = mcf; |
|
111 |
|
|
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 |
|
|
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 |
|
|
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 ⊃ |
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 |
|
|
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 |
|
|
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) |
... | ... |
@@ -199,136 +193,198 @@ |
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 |
|
210 |
PT result, bool optimal, typename CM::Value total, |
|
216 | 211 |
const std::string &test_id = "", |
217 |
|
|
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(); |
266 | 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; |
|
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, |
|
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( |
|
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, |
|
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, |
|
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, |
|
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). |
|
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 |
} |
... | ... |
@@ -9,27 +9,25 @@ |
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 |
... | ... |
@@ -110,26 +110,26 @@ |
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) |
0 comments (0 inline)