1.1 --- a/CMakeLists.txt Wed Apr 29 16:15:29 2009 +0100
1.2 +++ b/CMakeLists.txt Wed Apr 29 19:22:14 2009 +0100
1.3 @@ -17,8 +17,6 @@
1.4 FIND_PACKAGE(CPLEX)
1.5 FIND_PACKAGE(COIN)
1.6
1.7 -ADD_DEFINITIONS(-DHAVE_CONFIG_H)
1.8 -
1.9 IF(MSVC)
1.10 SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4250 /wd4355 /wd4800 /wd4996")
1.11 # Suppressed warnings:
1.12 @@ -28,8 +26,6 @@
1.13 # C4996: 'function': was declared deprecated
1.14 ENDIF(MSVC)
1.15
1.16 -ADD_DEFINITIONS(-DHAVE_CONFIG_H)
1.17 -
1.18 INCLUDE(CheckTypeSize)
1.19 CHECK_TYPE_SIZE("long long" LEMON_LONG_LONG)
1.20
2.1 --- a/Makefile.am Wed Apr 29 16:15:29 2009 +0100
2.2 +++ b/Makefile.am Wed Apr 29 19:22:14 2009 +0100
2.3 @@ -11,11 +11,12 @@
2.4 m4/lx_check_cplex.m4 \
2.5 m4/lx_check_glpk.m4 \
2.6 m4/lx_check_soplex.m4 \
2.7 - m4/lx_check_clp.m4 \
2.8 - m4/lx_check_cbc.m4 \
2.9 + m4/lx_check_coin.m4 \
2.10 CMakeLists.txt \
2.11 cmake/FindGhostscript.cmake \
2.12 + cmake/FindCPLEX.cmake \
2.13 cmake/FindGLPK.cmake \
2.14 + cmake/FindCOIN.cmake \
2.15 cmake/version.cmake.in \
2.16 cmake/version.cmake \
2.17 cmake/nsis/lemon.ico \
3.1 --- a/cmake/FindCOIN.cmake Wed Apr 29 16:15:29 2009 +0100
3.2 +++ b/cmake/FindCOIN.cmake Wed Apr 29 19:22:14 2009 +0100
3.3 @@ -1,28 +1,48 @@
3.4 SET(COIN_ROOT_DIR "" CACHE PATH "COIN root directory")
3.5
3.6 FIND_PATH(COIN_INCLUDE_DIR coin/CoinUtilsConfig.h
3.7 - PATHS ${COIN_ROOT_DIR}/include)
3.8 -
3.9 -FIND_LIBRARY(COIN_CBC_LIBRARY libCbc
3.10 - PATHS ${COIN_ROOT_DIR}/lib)
3.11 -FIND_LIBRARY(COIN_CBC_SOLVER_LIBRARY libCbcSolver
3.12 - PATHS ${COIN_ROOT_DIR}/lib)
3.13 -FIND_LIBRARY(COIN_CGL_LIBRARY libCgl
3.14 - PATHS ${COIN_ROOT_DIR}/lib)
3.15 -FIND_LIBRARY(COIN_CLP_LIBRARY libClp
3.16 - PATHS ${COIN_ROOT_DIR}/lib)
3.17 -FIND_LIBRARY(COIN_COIN_UTILS_LIBRARY libCoinUtils
3.18 - PATHS ${COIN_ROOT_DIR}/lib)
3.19 -FIND_LIBRARY(COIN_OSI_LIBRARY libOsi
3.20 - PATHS ${COIN_ROOT_DIR}/lib)
3.21 -FIND_LIBRARY(COIN_OSI_CBC_LIBRARY libOsiCbc
3.22 - PATHS ${COIN_ROOT_DIR}/lib)
3.23 -FIND_LIBRARY(COIN_OSI_CLP_LIBRARY libOsiClp
3.24 - PATHS ${COIN_ROOT_DIR}/lib)
3.25 -FIND_LIBRARY(COIN_OSI_VOL_LIBRARY libOsiVol
3.26 - PATHS ${COIN_ROOT_DIR}/lib)
3.27 -FIND_LIBRARY(COIN_VOL_LIBRARY libVol
3.28 - PATHS ${COIN_ROOT_DIR}/lib)
3.29 + HINTS ${COIN_ROOT_DIR}/include
3.30 +)
3.31 +FIND_LIBRARY(COIN_CBC_LIBRARY
3.32 + NAMES Cbc libCbc
3.33 + HINTS ${COIN_ROOT_DIR}/lib
3.34 +)
3.35 +FIND_LIBRARY(COIN_CBC_SOLVER_LIBRARY
3.36 + NAMES CbcSolver libCbcSolver
3.37 + HINTS ${COIN_ROOT_DIR}/lib
3.38 +)
3.39 +FIND_LIBRARY(COIN_CGL_LIBRARY
3.40 + NAMES Cgl libCgl
3.41 + HINTS ${COIN_ROOT_DIR}/lib
3.42 +)
3.43 +FIND_LIBRARY(COIN_CLP_LIBRARY
3.44 + NAMES Clp libClp
3.45 + HINTS ${COIN_ROOT_DIR}/lib
3.46 +)
3.47 +FIND_LIBRARY(COIN_COIN_UTILS_LIBRARY
3.48 + NAMES CoinUtils libCoinUtils
3.49 + HINTS ${COIN_ROOT_DIR}/lib
3.50 +)
3.51 +FIND_LIBRARY(COIN_OSI_LIBRARY
3.52 + NAMES Osi libOsi
3.53 + HINTS ${COIN_ROOT_DIR}/lib
3.54 +)
3.55 +FIND_LIBRARY(COIN_OSI_CBC_LIBRARY
3.56 + NAMES OsiCbc libOsiCbc
3.57 + HINTS ${COIN_ROOT_DIR}/lib
3.58 +)
3.59 +FIND_LIBRARY(COIN_OSI_CLP_LIBRARY
3.60 + NAMES OsiClp libOsiClp
3.61 + HINTS ${COIN_ROOT_DIR}/lib
3.62 +)
3.63 +FIND_LIBRARY(COIN_OSI_VOL_LIBRARY
3.64 + NAMES OsiVol libOsiVol
3.65 + HINTS ${COIN_ROOT_DIR}/lib
3.66 +)
3.67 +FIND_LIBRARY(COIN_VOL_LIBRARY
3.68 + NAMES Vol libVol
3.69 + HINTS ${COIN_ROOT_DIR}/lib
3.70 +)
3.71
3.72 INCLUDE(FindPackageHandleStandardArgs)
3.73 FIND_PACKAGE_HANDLE_STANDARD_ARGS(COIN DEFAULT_MSG
4.1 --- a/cmake/FindCPLEX.cmake Wed Apr 29 16:15:29 2009 +0100
4.2 +++ b/cmake/FindCPLEX.cmake Wed Apr 29 19:22:14 2009 +0100
4.3 @@ -1,21 +1,32 @@
4.4 +SET(CPLEX_ROOT_DIR "" CACHE PATH "CPLEX root directory")
4.5 +
4.6 FIND_PATH(CPLEX_INCLUDE_DIR
4.7 ilcplex/cplex.h
4.8 - PATHS "C:/ILOG/CPLEX91/include")
4.9 -
4.10 + PATHS "C:/ILOG/CPLEX91/include"
4.11 + PATHS "/opt/ilog/cplex91/include"
4.12 + HINTS ${CPLEX_ROOT_DIR}/include
4.13 +)
4.14 FIND_LIBRARY(CPLEX_LIBRARY
4.15 - NAMES cplex91
4.16 - PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda")
4.17 + cplex91
4.18 + PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda"
4.19 + PATHS "/opt/ilog/cplex91/bin"
4.20 + HINTS ${CPLEX_ROOT_DIR}/bin
4.21 +)
4.22
4.23 INCLUDE(FindPackageHandleStandardArgs)
4.24 FIND_PACKAGE_HANDLE_STANDARD_ARGS(CPLEX DEFAULT_MSG CPLEX_LIBRARY CPLEX_INCLUDE_DIR)
4.25
4.26 FIND_PATH(CPLEX_BIN_DIR
4.27 cplex91.dll
4.28 - PATHS "C:/ILOG/CPLEX91/bin/x86_win32")
4.29 + PATHS "C:/ILOG/CPLEX91/bin/x86_win32"
4.30 +)
4.31
4.32 IF(CPLEX_FOUND)
4.33 SET(CPLEX_INCLUDE_DIRS ${CPLEX_INCLUDE_DIR})
4.34 SET(CPLEX_LIBRARIES ${CPLEX_LIBRARY})
4.35 + IF(CMAKE_SYSTEM_NAME STREQUAL "Linux")
4.36 + SET(CPLEX_LIBRARIES "${CPLEX_LIBRARIES};m;pthread")
4.37 + ENDIF(CMAKE_SYSTEM_NAME STREQUAL "Linux")
4.38 ENDIF(CPLEX_FOUND)
4.39
4.40 MARK_AS_ADVANCED(CPLEX_LIBRARY CPLEX_INCLUDE_DIR CPLEX_BIN_DIR)
5.1 --- a/cmake/FindGLPK.cmake Wed Apr 29 16:15:29 2009 +0100
5.2 +++ b/cmake/FindGLPK.cmake Wed Apr 29 19:22:14 2009 +0100
5.3 @@ -1,16 +1,50 @@
5.4 +SET(GLPK_ROOT_DIR "" CACHE PATH "GLPK root directory")
5.5 +
5.6 SET(GLPK_REGKEY "[HKEY_LOCAL_MACHINE\\SOFTWARE\\GnuWin32\\Glpk;InstallPath]")
5.7 GET_FILENAME_COMPONENT(GLPK_ROOT_PATH ${GLPK_REGKEY} ABSOLUTE)
5.8
5.9 FIND_PATH(GLPK_INCLUDE_DIR
5.10 glpk.h
5.11 - PATHS ${GLPK_REGKEY}/include)
5.12 + PATHS ${GLPK_REGKEY}/include
5.13 + HINTS ${GLPK_ROOT_DIR}/include
5.14 +)
5.15 +FIND_LIBRARY(GLPK_LIBRARY
5.16 + glpk
5.17 + PATHS ${GLPK_REGKEY}/lib
5.18 + HINTS ${GLPK_ROOT_DIR}/lib
5.19 +)
5.20
5.21 -FIND_LIBRARY(GLPK_LIBRARY
5.22 - NAMES glpk
5.23 - PATHS ${GLPK_REGKEY}/lib)
5.24 +IF(GLPK_INCLUDE_DIR AND GLPK_LIBRARY)
5.25 + FILE(READ ${GLPK_INCLUDE_DIR}/glpk.h GLPK_GLPK_H)
5.26 +
5.27 + STRING(REGEX MATCH "define[ ]+GLP_MAJOR_VERSION[ ]+[0-9]+" GLPK_MAJOR_VERSION_LINE "${GLPK_GLPK_H}")
5.28 + STRING(REGEX REPLACE "define[ ]+GLP_MAJOR_VERSION[ ]+([0-9]+)" "\\1" GLPK_VERSION_MAJOR "${GLPK_MAJOR_VERSION_LINE}")
5.29 +
5.30 + STRING(REGEX MATCH "define[ ]+GLP_MINOR_VERSION[ ]+[0-9]+" GLPK_MINOR_VERSION_LINE "${GLPK_GLPK_H}")
5.31 + STRING(REGEX REPLACE "define[ ]+GLP_MINOR_VERSION[ ]+([0-9]+)" "\\1" GLPK_VERSION_MINOR "${GLPK_MINOR_VERSION_LINE}")
5.32 +
5.33 + SET(GLPK_VERSION_STRING "${GLPK_VERSION_MAJOR}.${GLPK_VERSION_MINOR}")
5.34 +
5.35 + IF(GLPK_FIND_VERSION)
5.36 + IF(GLPK_FIND_VERSION_COUNT GREATER 2)
5.37 + MESSAGE(SEND_ERROR "unexpected version string")
5.38 + ENDIF(GLPK_FIND_VERSION_COUNT GREATER 2)
5.39 +
5.40 + MATH(EXPR GLPK_REQUESTED_VERSION "${GLPK_FIND_VERSION_MAJOR}*100 + ${GLPK_FIND_VERSION_MINOR}")
5.41 + MATH(EXPR GLPK_FOUND_VERSION "${GLPK_VERSION_MAJOR}*100 + ${GLPK_VERSION_MINOR}")
5.42 +
5.43 + IF(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION)
5.44 + SET(GLPK_PROPER_VERSION_FOUND FALSE)
5.45 + ELSE(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION)
5.46 + SET(GLPK_PROPER_VERSION_FOUND TRUE)
5.47 + ENDIF(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION)
5.48 + ELSE(GLPK_FIND_VERSION)
5.49 + SET(GLPK_PROPER_VERSION_FOUND TRUE)
5.50 + ENDIF(GLPK_FIND_VERSION)
5.51 +ENDIF(GLPK_INCLUDE_DIR AND GLPK_LIBRARY)
5.52
5.53 INCLUDE(FindPackageHandleStandardArgs)
5.54 -FIND_PACKAGE_HANDLE_STANDARD_ARGS(GLPK DEFAULT_MSG GLPK_LIBRARY GLPK_INCLUDE_DIR)
5.55 +FIND_PACKAGE_HANDLE_STANDARD_ARGS(GLPK DEFAULT_MSG GLPK_LIBRARY GLPK_INCLUDE_DIR GLPK_PROPER_VERSION_FOUND)
5.56
5.57 IF(GLPK_FOUND)
5.58 SET(GLPK_INCLUDE_DIRS ${GLPK_INCLUDE_DIR})
6.1 --- a/doc/groups.dox Wed Apr 29 16:15:29 2009 +0100
6.2 +++ b/doc/groups.dox Wed Apr 29 19:22:14 2009 +0100
6.3 @@ -352,17 +352,17 @@
6.4 minimum total cost from a set of supply nodes to a set of demand nodes
6.5 in a network with capacity constraints (lower and upper bounds)
6.6 and arc costs.
6.7 -Formally, let \f$G=(V,A)\f$ be a digraph,
6.8 -\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
6.9 +Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
6.10 +\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
6.11 upper bounds for the flow values on the arcs, for which
6.12 -\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
6.13 -\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
6.14 -on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
6.15 +\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
6.16 +\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
6.17 +on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
6.18 signed supply values of the nodes.
6.19 If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
6.20 supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
6.21 \f$-sup(u)\f$ demand.
6.22 -A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
6.23 +A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
6.24 of the following optimization problem.
6.25
6.26 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
6.27 @@ -404,7 +404,7 @@
6.28
6.29 The dual solution of the minimum cost flow problem is represented by node
6.30 potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
6.31 -An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
6.32 +An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
6.33 is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
6.34 node potentials the following \e complementary \e slackness optimality
6.35 conditions hold.
6.36 @@ -413,15 +413,15 @@
6.37 - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
6.38 - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
6.39 - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
6.40 - - For all \f$u\in V\f$:
6.41 + - For all \f$u\in V\f$ nodes:
6.42 - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
6.43 then \f$\pi(u)=0\f$.
6.44
6.45 Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
6.46 -\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
6.47 +\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
6.48 \f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
6.49
6.50 -All algorithms provide dual solution (node potentials) as well
6.51 +All algorithms provide dual solution (node potentials) as well,
6.52 if an optimal flow is found.
6.53
6.54 LEMON contains several algorithms for solving minimum cost flow problems.
7.1 --- a/lemon/Makefile.am Wed Apr 29 16:15:29 2009 +0100
7.2 +++ b/lemon/Makefile.am Wed Apr 29 19:22:14 2009 +0100
7.3 @@ -15,7 +15,8 @@
7.4 lemon/random.cc \
7.5 lemon/bits/windows.cc
7.6
7.7 -
7.8 +nodist_lemon_HEADERS = lemon/config.h
7.9 +
7.10 lemon_libemon_la_CXXFLAGS = \
7.11 $(AM_CXXFLAGS) \
7.12 $(GLPK_CFLAGS) \
7.13 @@ -57,11 +58,11 @@
7.14 lemon/assert.h \
7.15 lemon/bfs.h \
7.16 lemon/bin_heap.h \
7.17 + lemon/cbc.h \
7.18 lemon/circulation.h \
7.19 lemon/clp.h \
7.20 lemon/color.h \
7.21 lemon/concept_check.h \
7.22 - lemon/config.h \
7.23 lemon/connectivity.h \
7.24 lemon/counter.h \
7.25 lemon/core.h \
8.1 --- a/lemon/circulation.h Wed Apr 29 16:15:29 2009 +0100
8.2 +++ b/lemon/circulation.h Wed Apr 29 19:22:14 2009 +0100
8.3 @@ -64,15 +64,15 @@
8.4 /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
8.5 typedef SM SupplyMap;
8.6
8.7 - /// \brief The type of the flow values.
8.8 - typedef typename SupplyMap::Value Flow;
8.9 + /// \brief The type of the flow and supply values.
8.10 + typedef typename SupplyMap::Value Value;
8.11
8.12 /// \brief The type of the map that stores the flow values.
8.13 ///
8.14 /// The type of the map that stores the flow values.
8.15 /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
8.16 /// concept.
8.17 - typedef typename Digraph::template ArcMap<Flow> FlowMap;
8.18 + typedef typename Digraph::template ArcMap<Value> FlowMap;
8.19
8.20 /// \brief Instantiates a FlowMap.
8.21 ///
8.22 @@ -104,7 +104,7 @@
8.23 /// \brief The tolerance used by the algorithm
8.24 ///
8.25 /// The tolerance used by the algorithm to handle inexact computation.
8.26 - typedef lemon::Tolerance<Flow> Tolerance;
8.27 + typedef lemon::Tolerance<Value> Tolerance;
8.28
8.29 };
8.30
8.31 @@ -187,8 +187,8 @@
8.32 typedef TR Traits;
8.33 ///The type of the digraph the algorithm runs on.
8.34 typedef typename Traits::Digraph Digraph;
8.35 - ///The type of the flow values.
8.36 - typedef typename Traits::Flow Flow;
8.37 + ///The type of the flow and supply values.
8.38 + typedef typename Traits::Value Value;
8.39
8.40 ///The type of the lower bound map.
8.41 typedef typename Traits::LowerMap LowerMap;
8.42 @@ -221,7 +221,7 @@
8.43 Elevator* _level;
8.44 bool _local_level;
8.45
8.46 - typedef typename Digraph::template NodeMap<Flow> ExcessMap;
8.47 + typedef typename Digraph::template NodeMap<Value> ExcessMap;
8.48 ExcessMap* _excess;
8.49
8.50 Tolerance _tol;
8.51 @@ -530,7 +530,7 @@
8.52 (*_excess)[_g.target(e)] += (*_lo)[e];
8.53 (*_excess)[_g.source(e)] -= (*_lo)[e];
8.54 } else {
8.55 - Flow fc = -(*_excess)[_g.target(e)];
8.56 + Value fc = -(*_excess)[_g.target(e)];
8.57 _flow->set(e, fc);
8.58 (*_excess)[_g.target(e)] = 0;
8.59 (*_excess)[_g.source(e)] -= fc;
8.60 @@ -563,11 +563,11 @@
8.61 while((act=_level->highestActive())!=INVALID) {
8.62 int actlevel=(*_level)[act];
8.63 int mlevel=_node_num;
8.64 - Flow exc=(*_excess)[act];
8.65 + Value exc=(*_excess)[act];
8.66
8.67 for(OutArcIt e(_g,act);e!=INVALID; ++e) {
8.68 Node v = _g.target(e);
8.69 - Flow fc=(*_up)[e]-(*_flow)[e];
8.70 + Value fc=(*_up)[e]-(*_flow)[e];
8.71 if(!_tol.positive(fc)) continue;
8.72 if((*_level)[v]<actlevel) {
8.73 if(!_tol.less(fc, exc)) {
8.74 @@ -591,7 +591,7 @@
8.75 }
8.76 for(InArcIt e(_g,act);e!=INVALID; ++e) {
8.77 Node v = _g.source(e);
8.78 - Flow fc=(*_flow)[e]-(*_lo)[e];
8.79 + Value fc=(*_flow)[e]-(*_lo)[e];
8.80 if(!_tol.positive(fc)) continue;
8.81 if((*_level)[v]<actlevel) {
8.82 if(!_tol.less(fc, exc)) {
8.83 @@ -661,13 +661,13 @@
8.84
8.85 ///@{
8.86
8.87 - /// \brief Returns the flow on the given arc.
8.88 + /// \brief Returns the flow value on the given arc.
8.89 ///
8.90 - /// Returns the flow on the given arc.
8.91 + /// Returns the flow value on the given arc.
8.92 ///
8.93 /// \pre Either \ref run() or \ref init() must be called before
8.94 /// using this function.
8.95 - Flow flow(const Arc& arc) const {
8.96 + Value flow(const Arc& arc) const {
8.97 return (*_flow)[arc];
8.98 }
8.99
8.100 @@ -750,7 +750,7 @@
8.101 if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
8.102 for(NodeIt n(_g);n!=INVALID;++n)
8.103 {
8.104 - Flow dif=-(*_supply)[n];
8.105 + Value dif=-(*_supply)[n];
8.106 for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
8.107 for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
8.108 if(_tol.negative(dif)) return false;
8.109 @@ -765,10 +765,10 @@
8.110 ///\sa barrierMap()
8.111 bool checkBarrier() const
8.112 {
8.113 - Flow delta=0;
8.114 - Flow inf_cap = std::numeric_limits<Flow>::has_infinity ?
8.115 - std::numeric_limits<Flow>::infinity() :
8.116 - std::numeric_limits<Flow>::max();
8.117 + Value delta=0;
8.118 + Value inf_cap = std::numeric_limits<Value>::has_infinity ?
8.119 + std::numeric_limits<Value>::infinity() :
8.120 + std::numeric_limits<Value>::max();
8.121 for(NodeIt n(_g);n!=INVALID;++n)
8.122 if(barrier(n))
8.123 delta-=(*_supply)[n];
9.1 --- a/lemon/core.h Wed Apr 29 16:15:29 2009 +0100
9.2 +++ b/lemon/core.h Wed Apr 29 19:22:14 2009 +0100
9.3 @@ -22,7 +22,7 @@
9.4 #include <vector>
9.5 #include <algorithm>
9.6
9.7 -#include <lemon/core.h>
9.8 +#include <lemon/config.h>
9.9 #include <lemon/bits/enable_if.h>
9.10 #include <lemon/bits/traits.h>
9.11 #include <lemon/assert.h>
10.1 --- a/lemon/network_simplex.h Wed Apr 29 16:15:29 2009 +0100
10.2 +++ b/lemon/network_simplex.h Wed Apr 29 19:22:14 2009 +0100
10.3 @@ -30,9 +30,6 @@
10.4
10.5 #include <lemon/core.h>
10.6 #include <lemon/math.h>
10.7 -#include <lemon/maps.h>
10.8 -#include <lemon/circulation.h>
10.9 -#include <lemon/adaptors.h>
10.10
10.11 namespace lemon {
10.12
10.13 @@ -50,14 +47,19 @@
10.14 ///
10.15 /// In general this class is the fastest implementation available
10.16 /// in LEMON for the minimum cost flow problem.
10.17 - /// Moreover it supports both direction of the supply/demand inequality
10.18 - /// constraints. For more information see \ref ProblemType.
10.19 + /// Moreover it supports both directions of the supply/demand inequality
10.20 + /// constraints. For more information see \ref SupplyType.
10.21 + ///
10.22 + /// Most of the parameters of the problem (except for the digraph)
10.23 + /// can be given using separate functions, and the algorithm can be
10.24 + /// executed using the \ref run() function. If some parameters are not
10.25 + /// specified, then default values will be used.
10.26 ///
10.27 /// \tparam GR The digraph type the algorithm runs on.
10.28 - /// \tparam F The value type used for flow amounts, capacity bounds
10.29 + /// \tparam V The value type used for flow amounts, capacity bounds
10.30 /// and supply values in the algorithm. By default it is \c int.
10.31 /// \tparam C The value type used for costs and potentials in the
10.32 - /// algorithm. By default it is the same as \c F.
10.33 + /// algorithm. By default it is the same as \c V.
10.34 ///
10.35 /// \warning Both value types must be signed and all input data must
10.36 /// be integer.
10.37 @@ -65,34 +67,92 @@
10.38 /// \note %NetworkSimplex provides five different pivot rule
10.39 /// implementations, from which the most efficient one is used
10.40 /// by default. For more information see \ref PivotRule.
10.41 - template <typename GR, typename F = int, typename C = F>
10.42 + template <typename GR, typename V = int, typename C = V>
10.43 class NetworkSimplex
10.44 {
10.45 public:
10.46
10.47 - /// The flow type of the algorithm
10.48 - typedef F Flow;
10.49 - /// The cost type of the algorithm
10.50 + /// The type of the flow amounts, capacity bounds and supply values
10.51 + typedef V Value;
10.52 + /// The type of the arc costs
10.53 typedef C Cost;
10.54 -#ifdef DOXYGEN
10.55 - /// The type of the flow map
10.56 - typedef GR::ArcMap<Flow> FlowMap;
10.57 - /// The type of the potential map
10.58 - typedef GR::NodeMap<Cost> PotentialMap;
10.59 -#else
10.60 - /// The type of the flow map
10.61 - typedef typename GR::template ArcMap<Flow> FlowMap;
10.62 - /// The type of the potential map
10.63 - typedef typename GR::template NodeMap<Cost> PotentialMap;
10.64 -#endif
10.65
10.66 public:
10.67
10.68 - /// \brief Enum type for selecting the pivot rule.
10.69 + /// \brief Problem type constants for the \c run() function.
10.70 ///
10.71 - /// Enum type for selecting the pivot rule for the \ref run()
10.72 + /// Enum type containing the problem type constants that can be
10.73 + /// returned by the \ref run() function of the algorithm.
10.74 + enum ProblemType {
10.75 + /// The problem has no feasible solution (flow).
10.76 + INFEASIBLE,
10.77 + /// The problem has optimal solution (i.e. it is feasible and
10.78 + /// bounded), and the algorithm has found optimal flow and node
10.79 + /// potentials (primal and dual solutions).
10.80 + OPTIMAL,
10.81 + /// The objective function of the problem is unbounded, i.e.
10.82 + /// there is a directed cycle having negative total cost and
10.83 + /// infinite upper bound.
10.84 + UNBOUNDED
10.85 + };
10.86 +
10.87 + /// \brief Constants for selecting the type of the supply constraints.
10.88 + ///
10.89 + /// Enum type containing constants for selecting the supply type,
10.90 + /// i.e. the direction of the inequalities in the supply/demand
10.91 + /// constraints of the \ref min_cost_flow "minimum cost flow problem".
10.92 + ///
10.93 + /// The default supply type is \c GEQ, since this form is supported
10.94 + /// by other minimum cost flow algorithms and the \ref Circulation
10.95 + /// algorithm, as well.
10.96 + /// The \c LEQ problem type can be selected using the \ref supplyType()
10.97 /// function.
10.98 ///
10.99 + /// Note that the equality form is a special case of both supply types.
10.100 + enum SupplyType {
10.101 +
10.102 + /// This option means that there are <em>"greater or equal"</em>
10.103 + /// supply/demand constraints in the definition, i.e. the exact
10.104 + /// formulation of the problem is the following.
10.105 + /**
10.106 + \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
10.107 + \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
10.108 + sup(u) \quad \forall u\in V \f]
10.109 + \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
10.110 + */
10.111 + /// It means that the total demand must be greater or equal to the
10.112 + /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
10.113 + /// negative) and all the supplies have to be carried out from
10.114 + /// the supply nodes, but there could be demands that are not
10.115 + /// satisfied.
10.116 + GEQ,
10.117 + /// It is just an alias for the \c GEQ option.
10.118 + CARRY_SUPPLIES = GEQ,
10.119 +
10.120 + /// This option means that there are <em>"less or equal"</em>
10.121 + /// supply/demand constraints in the definition, i.e. the exact
10.122 + /// formulation of the problem is the following.
10.123 + /**
10.124 + \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
10.125 + \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
10.126 + sup(u) \quad \forall u\in V \f]
10.127 + \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
10.128 + */
10.129 + /// It means that the total demand must be less or equal to the
10.130 + /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
10.131 + /// positive) and all the demands have to be satisfied, but there
10.132 + /// could be supplies that are not carried out from the supply
10.133 + /// nodes.
10.134 + LEQ,
10.135 + /// It is just an alias for the \c LEQ option.
10.136 + SATISFY_DEMANDS = LEQ
10.137 + };
10.138 +
10.139 + /// \brief Constants for selecting the pivot rule.
10.140 + ///
10.141 + /// Enum type containing constants for selecting the pivot rule for
10.142 + /// the \ref run() function.
10.143 + ///
10.144 /// \ref NetworkSimplex provides five different pivot rule
10.145 /// implementations that significantly affect the running time
10.146 /// of the algorithm.
10.147 @@ -131,71 +191,15 @@
10.148 ALTERING_LIST
10.149 };
10.150
10.151 - /// \brief Enum type for selecting the problem type.
10.152 - ///
10.153 - /// Enum type for selecting the problem type, i.e. the direction of
10.154 - /// the inequalities in the supply/demand constraints of the
10.155 - /// \ref min_cost_flow "minimum cost flow problem".
10.156 - ///
10.157 - /// The default problem type is \c GEQ, since this form is supported
10.158 - /// by other minimum cost flow algorithms and the \ref Circulation
10.159 - /// algorithm as well.
10.160 - /// The \c LEQ problem type can be selected using the \ref problemType()
10.161 - /// function.
10.162 - ///
10.163 - /// Note that the equality form is a special case of both problem type.
10.164 - enum ProblemType {
10.165 -
10.166 - /// This option means that there are "<em>greater or equal</em>"
10.167 - /// constraints in the defintion, i.e. the exact formulation of the
10.168 - /// problem is the following.
10.169 - /**
10.170 - \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
10.171 - \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
10.172 - sup(u) \quad \forall u\in V \f]
10.173 - \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
10.174 - */
10.175 - /// It means that the total demand must be greater or equal to the
10.176 - /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
10.177 - /// negative) and all the supplies have to be carried out from
10.178 - /// the supply nodes, but there could be demands that are not
10.179 - /// satisfied.
10.180 - GEQ,
10.181 - /// It is just an alias for the \c GEQ option.
10.182 - CARRY_SUPPLIES = GEQ,
10.183 -
10.184 - /// This option means that there are "<em>less or equal</em>"
10.185 - /// constraints in the defintion, i.e. the exact formulation of the
10.186 - /// problem is the following.
10.187 - /**
10.188 - \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
10.189 - \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
10.190 - sup(u) \quad \forall u\in V \f]
10.191 - \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
10.192 - */
10.193 - /// It means that the total demand must be less or equal to the
10.194 - /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
10.195 - /// positive) and all the demands have to be satisfied, but there
10.196 - /// could be supplies that are not carried out from the supply
10.197 - /// nodes.
10.198 - LEQ,
10.199 - /// It is just an alias for the \c LEQ option.
10.200 - SATISFY_DEMANDS = LEQ
10.201 - };
10.202 -
10.203 private:
10.204
10.205 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
10.206
10.207 - typedef typename GR::template ArcMap<Flow> FlowArcMap;
10.208 - typedef typename GR::template ArcMap<Cost> CostArcMap;
10.209 - typedef typename GR::template NodeMap<Flow> FlowNodeMap;
10.210 -
10.211 typedef std::vector<Arc> ArcVector;
10.212 typedef std::vector<Node> NodeVector;
10.213 typedef std::vector<int> IntVector;
10.214 typedef std::vector<bool> BoolVector;
10.215 - typedef std::vector<Flow> FlowVector;
10.216 + typedef std::vector<Value> ValueVector;
10.217 typedef std::vector<Cost> CostVector;
10.218
10.219 // State constants for arcs
10.220 @@ -213,32 +217,23 @@
10.221 int _arc_num;
10.222
10.223 // Parameters of the problem
10.224 - FlowArcMap *_plower;
10.225 - FlowArcMap *_pupper;
10.226 - CostArcMap *_pcost;
10.227 - FlowNodeMap *_psupply;
10.228 - bool _pstsup;
10.229 - Node _psource, _ptarget;
10.230 - Flow _pstflow;
10.231 - ProblemType _ptype;
10.232 -
10.233 - // Result maps
10.234 - FlowMap *_flow_map;
10.235 - PotentialMap *_potential_map;
10.236 - bool _local_flow;
10.237 - bool _local_potential;
10.238 + bool _have_lower;
10.239 + SupplyType _stype;
10.240 + Value _sum_supply;
10.241
10.242 // Data structures for storing the digraph
10.243 IntNodeMap _node_id;
10.244 - ArcVector _arc_ref;
10.245 + IntArcMap _arc_id;
10.246 IntVector _source;
10.247 IntVector _target;
10.248
10.249 // Node and arc data
10.250 - FlowVector _cap;
10.251 + ValueVector _lower;
10.252 + ValueVector _upper;
10.253 + ValueVector _cap;
10.254 CostVector _cost;
10.255 - FlowVector _supply;
10.256 - FlowVector _flow;
10.257 + ValueVector _supply;
10.258 + ValueVector _flow;
10.259 CostVector _pi;
10.260
10.261 // Data for storing the spanning tree structure
10.262 @@ -257,7 +252,16 @@
10.263 int in_arc, join, u_in, v_in, u_out, v_out;
10.264 int first, second, right, last;
10.265 int stem, par_stem, new_stem;
10.266 - Flow delta;
10.267 + Value delta;
10.268 +
10.269 + public:
10.270 +
10.271 + /// \brief Constant for infinite upper bounds (capacities).
10.272 + ///
10.273 + /// Constant for infinite upper bounds (capacities).
10.274 + /// It is \c std::numeric_limits<Value>::infinity() if available,
10.275 + /// \c std::numeric_limits<Value>::max() otherwise.
10.276 + const Value INF;
10.277
10.278 private:
10.279
10.280 @@ -659,25 +663,68 @@
10.281 ///
10.282 /// \param graph The digraph the algorithm runs on.
10.283 NetworkSimplex(const GR& graph) :
10.284 - _graph(graph),
10.285 - _plower(NULL), _pupper(NULL), _pcost(NULL),
10.286 - _psupply(NULL), _pstsup(false), _ptype(GEQ),
10.287 - _flow_map(NULL), _potential_map(NULL),
10.288 - _local_flow(false), _local_potential(false),
10.289 - _node_id(graph)
10.290 + _graph(graph), _node_id(graph), _arc_id(graph),
10.291 + INF(std::numeric_limits<Value>::has_infinity ?
10.292 + std::numeric_limits<Value>::infinity() :
10.293 + std::numeric_limits<Value>::max())
10.294 {
10.295 - LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
10.296 - std::numeric_limits<Flow>::is_signed,
10.297 - "The flow type of NetworkSimplex must be signed integer");
10.298 - LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
10.299 - std::numeric_limits<Cost>::is_signed,
10.300 - "The cost type of NetworkSimplex must be signed integer");
10.301 - }
10.302 + // Check the value types
10.303 + LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
10.304 + "The flow type of NetworkSimplex must be signed");
10.305 + LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
10.306 + "The cost type of NetworkSimplex must be signed");
10.307 +
10.308 + // Resize vectors
10.309 + _node_num = countNodes(_graph);
10.310 + _arc_num = countArcs(_graph);
10.311 + int all_node_num = _node_num + 1;
10.312 + int all_arc_num = _arc_num + _node_num;
10.313
10.314 - /// Destructor.
10.315 - ~NetworkSimplex() {
10.316 - if (_local_flow) delete _flow_map;
10.317 - if (_local_potential) delete _potential_map;
10.318 + _source.resize(all_arc_num);
10.319 + _target.resize(all_arc_num);
10.320 +
10.321 + _lower.resize(all_arc_num);
10.322 + _upper.resize(all_arc_num);
10.323 + _cap.resize(all_arc_num);
10.324 + _cost.resize(all_arc_num);
10.325 + _supply.resize(all_node_num);
10.326 + _flow.resize(all_arc_num);
10.327 + _pi.resize(all_node_num);
10.328 +
10.329 + _parent.resize(all_node_num);
10.330 + _pred.resize(all_node_num);
10.331 + _forward.resize(all_node_num);
10.332 + _thread.resize(all_node_num);
10.333 + _rev_thread.resize(all_node_num);
10.334 + _succ_num.resize(all_node_num);
10.335 + _last_succ.resize(all_node_num);
10.336 + _state.resize(all_arc_num);
10.337 +
10.338 + // Copy the graph (store the arcs in a mixed order)
10.339 + int i = 0;
10.340 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
10.341 + _node_id[n] = i;
10.342 + }
10.343 + int k = std::max(int(std::sqrt(double(_arc_num))), 10);
10.344 + i = 0;
10.345 + for (ArcIt a(_graph); a != INVALID; ++a) {
10.346 + _arc_id[a] = i;
10.347 + _source[i] = _node_id[_graph.source(a)];
10.348 + _target[i] = _node_id[_graph.target(a)];
10.349 + if ((i += k) >= _arc_num) i = (i % k) + 1;
10.350 + }
10.351 +
10.352 + // Initialize maps
10.353 + for (int i = 0; i != _node_num; ++i) {
10.354 + _supply[i] = 0;
10.355 + }
10.356 + for (int i = 0; i != _arc_num; ++i) {
10.357 + _lower[i] = 0;
10.358 + _upper[i] = INF;
10.359 + _cost[i] = 1;
10.360 + }
10.361 + _have_lower = false;
10.362 + _stype = GEQ;
10.363 }
10.364
10.365 /// \name Parameters
10.366 @@ -689,21 +736,19 @@
10.367 /// \brief Set the lower bounds on the arcs.
10.368 ///
10.369 /// This function sets the lower bounds on the arcs.
10.370 - /// If neither this function nor \ref boundMaps() is used before
10.371 - /// calling \ref run(), the lower bounds will be set to zero
10.372 - /// on all arcs.
10.373 + /// If it is not used before calling \ref run(), the lower bounds
10.374 + /// will be set to zero on all arcs.
10.375 ///
10.376 /// \param map An arc map storing the lower bounds.
10.377 - /// Its \c Value type must be convertible to the \c Flow type
10.378 + /// Its \c Value type must be convertible to the \c Value type
10.379 /// of the algorithm.
10.380 ///
10.381 /// \return <tt>(*this)</tt>
10.382 - template <typename LOWER>
10.383 - NetworkSimplex& lowerMap(const LOWER& map) {
10.384 - delete _plower;
10.385 - _plower = new FlowArcMap(_graph);
10.386 + template <typename LowerMap>
10.387 + NetworkSimplex& lowerMap(const LowerMap& map) {
10.388 + _have_lower = true;
10.389 for (ArcIt a(_graph); a != INVALID; ++a) {
10.390 - (*_plower)[a] = map[a];
10.391 + _lower[_arc_id[a]] = map[a];
10.392 }
10.393 return *this;
10.394 }
10.395 @@ -711,63 +756,23 @@
10.396 /// \brief Set the upper bounds (capacities) on the arcs.
10.397 ///
10.398 /// This function sets the upper bounds (capacities) on the arcs.
10.399 - /// If none of the functions \ref upperMap(), \ref capacityMap()
10.400 - /// and \ref boundMaps() is used before calling \ref run(),
10.401 - /// the upper bounds (capacities) will be set to
10.402 - /// \c std::numeric_limits<Flow>::max() on all arcs.
10.403 + /// If it is not used before calling \ref run(), the upper bounds
10.404 + /// will be set to \ref INF on all arcs (i.e. the flow value will be
10.405 + /// unbounded from above on each arc).
10.406 ///
10.407 /// \param map An arc map storing the upper bounds.
10.408 - /// Its \c Value type must be convertible to the \c Flow type
10.409 + /// Its \c Value type must be convertible to the \c Value type
10.410 /// of the algorithm.
10.411 ///
10.412 /// \return <tt>(*this)</tt>
10.413 - template<typename UPPER>
10.414 - NetworkSimplex& upperMap(const UPPER& map) {
10.415 - delete _pupper;
10.416 - _pupper = new FlowArcMap(_graph);
10.417 + template<typename UpperMap>
10.418 + NetworkSimplex& upperMap(const UpperMap& map) {
10.419 for (ArcIt a(_graph); a != INVALID; ++a) {
10.420 - (*_pupper)[a] = map[a];
10.421 + _upper[_arc_id[a]] = map[a];
10.422 }
10.423 return *this;
10.424 }
10.425
10.426 - /// \brief Set the upper bounds (capacities) on the arcs.
10.427 - ///
10.428 - /// This function sets the upper bounds (capacities) on the arcs.
10.429 - /// It is just an alias for \ref upperMap().
10.430 - ///
10.431 - /// \return <tt>(*this)</tt>
10.432 - template<typename CAP>
10.433 - NetworkSimplex& capacityMap(const CAP& map) {
10.434 - return upperMap(map);
10.435 - }
10.436 -
10.437 - /// \brief Set the lower and upper bounds on the arcs.
10.438 - ///
10.439 - /// This function sets the lower and upper bounds on the arcs.
10.440 - /// If neither this function nor \ref lowerMap() is used before
10.441 - /// calling \ref run(), the lower bounds will be set to zero
10.442 - /// on all arcs.
10.443 - /// If none of the functions \ref upperMap(), \ref capacityMap()
10.444 - /// and \ref boundMaps() is used before calling \ref run(),
10.445 - /// the upper bounds (capacities) will be set to
10.446 - /// \c std::numeric_limits<Flow>::max() on all arcs.
10.447 - ///
10.448 - /// \param lower An arc map storing the lower bounds.
10.449 - /// \param upper An arc map storing the upper bounds.
10.450 - ///
10.451 - /// The \c Value type of the maps must be convertible to the
10.452 - /// \c Flow type of the algorithm.
10.453 - ///
10.454 - /// \note This function is just a shortcut of calling \ref lowerMap()
10.455 - /// and \ref upperMap() separately.
10.456 - ///
10.457 - /// \return <tt>(*this)</tt>
10.458 - template <typename LOWER, typename UPPER>
10.459 - NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
10.460 - return lowerMap(lower).upperMap(upper);
10.461 - }
10.462 -
10.463 /// \brief Set the costs of the arcs.
10.464 ///
10.465 /// This function sets the costs of the arcs.
10.466 @@ -779,12 +784,10 @@
10.467 /// of the algorithm.
10.468 ///
10.469 /// \return <tt>(*this)</tt>
10.470 - template<typename COST>
10.471 - NetworkSimplex& costMap(const COST& map) {
10.472 - delete _pcost;
10.473 - _pcost = new CostArcMap(_graph);
10.474 + template<typename CostMap>
10.475 + NetworkSimplex& costMap(const CostMap& map) {
10.476 for (ArcIt a(_graph); a != INVALID; ++a) {
10.477 - (*_pcost)[a] = map[a];
10.478 + _cost[_arc_id[a]] = map[a];
10.479 }
10.480 return *this;
10.481 }
10.482 @@ -797,17 +800,14 @@
10.483 /// (It makes sense only if non-zero lower bounds are given.)
10.484 ///
10.485 /// \param map A node map storing the supply values.
10.486 - /// Its \c Value type must be convertible to the \c Flow type
10.487 + /// Its \c Value type must be convertible to the \c Value type
10.488 /// of the algorithm.
10.489 ///
10.490 /// \return <tt>(*this)</tt>
10.491 - template<typename SUP>
10.492 - NetworkSimplex& supplyMap(const SUP& map) {
10.493 - delete _psupply;
10.494 - _pstsup = false;
10.495 - _psupply = new FlowNodeMap(_graph);
10.496 + template<typename SupplyMap>
10.497 + NetworkSimplex& supplyMap(const SupplyMap& map) {
10.498 for (NodeIt n(_graph); n != INVALID; ++n) {
10.499 - (*_psupply)[n] = map[n];
10.500 + _supply[_node_id[n]] = map[n];
10.501 }
10.502 return *this;
10.503 }
10.504 @@ -820,71 +820,39 @@
10.505 /// calling \ref run(), the supply of each node will be set to zero.
10.506 /// (It makes sense only if non-zero lower bounds are given.)
10.507 ///
10.508 + /// Using this function has the same effect as using \ref supplyMap()
10.509 + /// with such a map in which \c k is assigned to \c s, \c -k is
10.510 + /// assigned to \c t and all other nodes have zero supply value.
10.511 + ///
10.512 /// \param s The source node.
10.513 /// \param t The target node.
10.514 /// \param k The required amount of flow from node \c s to node \c t
10.515 /// (i.e. the supply of \c s and the demand of \c t).
10.516 ///
10.517 /// \return <tt>(*this)</tt>
10.518 - NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
10.519 - delete _psupply;
10.520 - _psupply = NULL;
10.521 - _pstsup = true;
10.522 - _psource = s;
10.523 - _ptarget = t;
10.524 - _pstflow = k;
10.525 + NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
10.526 + for (int i = 0; i != _node_num; ++i) {
10.527 + _supply[i] = 0;
10.528 + }
10.529 + _supply[_node_id[s]] = k;
10.530 + _supply[_node_id[t]] = -k;
10.531 return *this;
10.532 }
10.533
10.534 - /// \brief Set the problem type.
10.535 + /// \brief Set the type of the supply constraints.
10.536 ///
10.537 - /// This function sets the problem type for the algorithm.
10.538 - /// If it is not used before calling \ref run(), the \ref GEQ problem
10.539 + /// This function sets the type of the supply/demand constraints.
10.540 + /// If it is not used before calling \ref run(), the \ref GEQ supply
10.541 /// type will be used.
10.542 ///
10.543 - /// For more information see \ref ProblemType.
10.544 + /// For more information see \ref SupplyType.
10.545 ///
10.546 /// \return <tt>(*this)</tt>
10.547 - NetworkSimplex& problemType(ProblemType problem_type) {
10.548 - _ptype = problem_type;
10.549 + NetworkSimplex& supplyType(SupplyType supply_type) {
10.550 + _stype = supply_type;
10.551 return *this;
10.552 }
10.553
10.554 - /// \brief Set the flow map.
10.555 - ///
10.556 - /// This function sets the flow map.
10.557 - /// If it is not used before calling \ref run(), an instance will
10.558 - /// be allocated automatically. The destructor deallocates this
10.559 - /// automatically allocated map, of course.
10.560 - ///
10.561 - /// \return <tt>(*this)</tt>
10.562 - NetworkSimplex& flowMap(FlowMap& map) {
10.563 - if (_local_flow) {
10.564 - delete _flow_map;
10.565 - _local_flow = false;
10.566 - }
10.567 - _flow_map = ↦
10.568 - return *this;
10.569 - }
10.570 -
10.571 - /// \brief Set the potential map.
10.572 - ///
10.573 - /// This function sets the potential map, which is used for storing
10.574 - /// the dual solution.
10.575 - /// If it is not used before calling \ref run(), an instance will
10.576 - /// be allocated automatically. The destructor deallocates this
10.577 - /// automatically allocated map, of course.
10.578 - ///
10.579 - /// \return <tt>(*this)</tt>
10.580 - NetworkSimplex& potentialMap(PotentialMap& map) {
10.581 - if (_local_potential) {
10.582 - delete _potential_map;
10.583 - _local_potential = false;
10.584 - }
10.585 - _potential_map = ↦
10.586 - return *this;
10.587 - }
10.588 -
10.589 /// @}
10.590
10.591 /// \name Execution Control
10.592 @@ -896,13 +864,12 @@
10.593 ///
10.594 /// This function runs the algorithm.
10.595 /// The paramters can be specified using functions \ref lowerMap(),
10.596 - /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
10.597 - /// \ref costMap(), \ref supplyMap(), \ref stSupply(),
10.598 - /// \ref problemType(), \ref flowMap() and \ref potentialMap().
10.599 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
10.600 + /// \ref supplyType().
10.601 /// For example,
10.602 /// \code
10.603 /// NetworkSimplex<ListDigraph> ns(graph);
10.604 - /// ns.boundMaps(lower, upper).costMap(cost)
10.605 + /// ns.lowerMap(lower).upperMap(upper).costMap(cost)
10.606 /// .supplyMap(sup).run();
10.607 /// \endcode
10.608 ///
10.609 @@ -910,33 +877,44 @@
10.610 /// that have been given are kept for the next call, unless
10.611 /// \ref reset() is called, thus only the modified parameters
10.612 /// have to be set again. See \ref reset() for examples.
10.613 + /// However the underlying digraph must not be modified after this
10.614 + /// class have been constructed, since it copies and extends the graph.
10.615 ///
10.616 /// \param pivot_rule The pivot rule that will be used during the
10.617 /// algorithm. For more information see \ref PivotRule.
10.618 ///
10.619 - /// \return \c true if a feasible flow can be found.
10.620 - bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
10.621 - return init() && start(pivot_rule);
10.622 + /// \return \c INFEASIBLE if no feasible flow exists,
10.623 + /// \n \c OPTIMAL if the problem has optimal solution
10.624 + /// (i.e. it is feasible and bounded), and the algorithm has found
10.625 + /// optimal flow and node potentials (primal and dual solutions),
10.626 + /// \n \c UNBOUNDED if the objective function of the problem is
10.627 + /// unbounded, i.e. there is a directed cycle having negative total
10.628 + /// cost and infinite upper bound.
10.629 + ///
10.630 + /// \see ProblemType, PivotRule
10.631 + ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
10.632 + if (!init()) return INFEASIBLE;
10.633 + return start(pivot_rule);
10.634 }
10.635
10.636 /// \brief Reset all the parameters that have been given before.
10.637 ///
10.638 /// This function resets all the paramaters that have been given
10.639 /// before using functions \ref lowerMap(), \ref upperMap(),
10.640 - /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
10.641 - /// \ref supplyMap(), \ref stSupply(), \ref problemType(),
10.642 - /// \ref flowMap() and \ref potentialMap().
10.643 + /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
10.644 ///
10.645 /// It is useful for multiple run() calls. If this function is not
10.646 /// used, all the parameters given before are kept for the next
10.647 /// \ref run() call.
10.648 + /// However the underlying digraph must not be modified after this
10.649 + /// class have been constructed, since it copies and extends the graph.
10.650 ///
10.651 /// For example,
10.652 /// \code
10.653 /// NetworkSimplex<ListDigraph> ns(graph);
10.654 ///
10.655 /// // First run
10.656 - /// ns.lowerMap(lower).capacityMap(cap).costMap(cost)
10.657 + /// ns.lowerMap(lower).upperMap(upper).costMap(cost)
10.658 /// .supplyMap(sup).run();
10.659 ///
10.660 /// // Run again with modified cost map (reset() is not called,
10.661 @@ -947,29 +925,22 @@
10.662 /// // Run again from scratch using reset()
10.663 /// // (the lower bounds will be set to zero on all arcs)
10.664 /// ns.reset();
10.665 - /// ns.capacityMap(cap).costMap(cost)
10.666 + /// ns.upperMap(capacity).costMap(cost)
10.667 /// .supplyMap(sup).run();
10.668 /// \endcode
10.669 ///
10.670 /// \return <tt>(*this)</tt>
10.671 NetworkSimplex& reset() {
10.672 - delete _plower;
10.673 - delete _pupper;
10.674 - delete _pcost;
10.675 - delete _psupply;
10.676 - _plower = NULL;
10.677 - _pupper = NULL;
10.678 - _pcost = NULL;
10.679 - _psupply = NULL;
10.680 - _pstsup = false;
10.681 - _ptype = GEQ;
10.682 - if (_local_flow) delete _flow_map;
10.683 - if (_local_potential) delete _potential_map;
10.684 - _flow_map = NULL;
10.685 - _potential_map = NULL;
10.686 - _local_flow = false;
10.687 - _local_potential = false;
10.688 -
10.689 + for (int i = 0; i != _node_num; ++i) {
10.690 + _supply[i] = 0;
10.691 + }
10.692 + for (int i = 0; i != _arc_num; ++i) {
10.693 + _lower[i] = 0;
10.694 + _upper[i] = INF;
10.695 + _cost[i] = 1;
10.696 + }
10.697 + _have_lower = false;
10.698 + _stype = GEQ;
10.699 return *this;
10.700 }
10.701
10.702 @@ -985,7 +956,7 @@
10.703 /// \brief Return the total cost of the found flow.
10.704 ///
10.705 /// This function returns the total cost of the found flow.
10.706 - /// The complexity of the function is O(e).
10.707 + /// Its complexity is O(e).
10.708 ///
10.709 /// \note The return type of the function can be specified as a
10.710 /// template parameter. For example,
10.711 @@ -997,15 +968,12 @@
10.712 /// function.
10.713 ///
10.714 /// \pre \ref run() must be called before using this function.
10.715 - template <typename Num>
10.716 - Num totalCost() const {
10.717 - Num c = 0;
10.718 - if (_pcost) {
10.719 - for (ArcIt e(_graph); e != INVALID; ++e)
10.720 - c += (*_flow_map)[e] * (*_pcost)[e];
10.721 - } else {
10.722 - for (ArcIt e(_graph); e != INVALID; ++e)
10.723 - c += (*_flow_map)[e];
10.724 + template <typename Number>
10.725 + Number totalCost() const {
10.726 + Number c = 0;
10.727 + for (ArcIt a(_graph); a != INVALID; ++a) {
10.728 + int i = _arc_id[a];
10.729 + c += Number(_flow[i]) * Number(_cost[i]);
10.730 }
10.731 return c;
10.732 }
10.733 @@ -1021,18 +989,22 @@
10.734 /// This function returns the flow on the given arc.
10.735 ///
10.736 /// \pre \ref run() must be called before using this function.
10.737 - Flow flow(const Arc& a) const {
10.738 - return (*_flow_map)[a];
10.739 + Value flow(const Arc& a) const {
10.740 + return _flow[_arc_id[a]];
10.741 }
10.742
10.743 - /// \brief Return a const reference to the flow map.
10.744 + /// \brief Return the flow map (the primal solution).
10.745 ///
10.746 - /// This function returns a const reference to an arc map storing
10.747 - /// the found flow.
10.748 + /// This function copies the flow value on each arc into the given
10.749 + /// map. The \c Value type of the algorithm must be convertible to
10.750 + /// the \c Value type of the map.
10.751 ///
10.752 /// \pre \ref run() must be called before using this function.
10.753 - const FlowMap& flowMap() const {
10.754 - return *_flow_map;
10.755 + template <typename FlowMap>
10.756 + void flowMap(FlowMap &map) const {
10.757 + for (ArcIt a(_graph); a != INVALID; ++a) {
10.758 + map.set(a, _flow[_arc_id[a]]);
10.759 + }
10.760 }
10.761
10.762 /// \brief Return the potential (dual value) of the given node.
10.763 @@ -1042,19 +1014,22 @@
10.764 ///
10.765 /// \pre \ref run() must be called before using this function.
10.766 Cost potential(const Node& n) const {
10.767 - return (*_potential_map)[n];
10.768 + return _pi[_node_id[n]];
10.769 }
10.770
10.771 - /// \brief Return a const reference to the potential map
10.772 - /// (the dual solution).
10.773 + /// \brief Return the potential map (the dual solution).
10.774 ///
10.775 - /// This function returns a const reference to a node map storing
10.776 - /// the found potentials, which form the dual solution of the
10.777 - /// \ref min_cost_flow "minimum cost flow" problem.
10.778 + /// This function copies the potential (dual value) of each node
10.779 + /// into the given map.
10.780 + /// The \c Cost type of the algorithm must be convertible to the
10.781 + /// \c Value type of the map.
10.782 ///
10.783 /// \pre \ref run() must be called before using this function.
10.784 - const PotentialMap& potentialMap() const {
10.785 - return *_potential_map;
10.786 + template <typename PotentialMap>
10.787 + void potentialMap(PotentialMap &map) const {
10.788 + for (NodeIt n(_graph); n != INVALID; ++n) {
10.789 + map.set(n, _pi[_node_id[n]]);
10.790 + }
10.791 }
10.792
10.793 /// @}
10.794 @@ -1063,245 +1038,82 @@
10.795
10.796 // Initialize internal data structures
10.797 bool init() {
10.798 - // Initialize result maps
10.799 - if (!_flow_map) {
10.800 - _flow_map = new FlowMap(_graph);
10.801 - _local_flow = true;
10.802 + if (_node_num == 0) return false;
10.803 +
10.804 + // Check the sum of supply values
10.805 + _sum_supply = 0;
10.806 + for (int i = 0; i != _node_num; ++i) {
10.807 + _sum_supply += _supply[i];
10.808 }
10.809 - if (!_potential_map) {
10.810 - _potential_map = new PotentialMap(_graph);
10.811 - _local_potential = true;
10.812 + if ( !((_stype == GEQ && _sum_supply <= 0) ||
10.813 + (_stype == LEQ && _sum_supply >= 0)) ) return false;
10.814 +
10.815 + // Remove non-zero lower bounds
10.816 + if (_have_lower) {
10.817 + for (int i = 0; i != _arc_num; ++i) {
10.818 + Value c = _lower[i];
10.819 + if (c >= 0) {
10.820 + _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
10.821 + } else {
10.822 + _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
10.823 + }
10.824 + _supply[_source[i]] -= c;
10.825 + _supply[_target[i]] += c;
10.826 + }
10.827 + } else {
10.828 + for (int i = 0; i != _arc_num; ++i) {
10.829 + _cap[i] = _upper[i];
10.830 + }
10.831 }
10.832
10.833 - // Initialize vectors
10.834 - _node_num = countNodes(_graph);
10.835 - _arc_num = countArcs(_graph);
10.836 - int all_node_num = _node_num + 1;
10.837 - int all_arc_num = _arc_num + _node_num;
10.838 - if (_node_num == 0) return false;
10.839 -
10.840 - _arc_ref.resize(_arc_num);
10.841 - _source.resize(all_arc_num);
10.842 - _target.resize(all_arc_num);
10.843 -
10.844 - _cap.resize(all_arc_num);
10.845 - _cost.resize(all_arc_num);
10.846 - _supply.resize(all_node_num);
10.847 - _flow.resize(all_arc_num);
10.848 - _pi.resize(all_node_num);
10.849 -
10.850 - _parent.resize(all_node_num);
10.851 - _pred.resize(all_node_num);
10.852 - _forward.resize(all_node_num);
10.853 - _thread.resize(all_node_num);
10.854 - _rev_thread.resize(all_node_num);
10.855 - _succ_num.resize(all_node_num);
10.856 - _last_succ.resize(all_node_num);
10.857 - _state.resize(all_arc_num);
10.858 -
10.859 - // Initialize node related data
10.860 - bool valid_supply = true;
10.861 - Flow sum_supply = 0;
10.862 - if (!_pstsup && !_psupply) {
10.863 - _pstsup = true;
10.864 - _psource = _ptarget = NodeIt(_graph);
10.865 - _pstflow = 0;
10.866 - }
10.867 - if (_psupply) {
10.868 - int i = 0;
10.869 - for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
10.870 - _node_id[n] = i;
10.871 - _supply[i] = (*_psupply)[n];
10.872 - sum_supply += _supply[i];
10.873 + // Initialize artifical cost
10.874 + Cost ART_COST;
10.875 + if (std::numeric_limits<Cost>::is_exact) {
10.876 + ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
10.877 + } else {
10.878 + ART_COST = std::numeric_limits<Cost>::min();
10.879 + for (int i = 0; i != _arc_num; ++i) {
10.880 + if (_cost[i] > ART_COST) ART_COST = _cost[i];
10.881 }
10.882 - valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
10.883 - (_ptype == LEQ && sum_supply >= 0);
10.884 - } else {
10.885 - int i = 0;
10.886 - for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
10.887 - _node_id[n] = i;
10.888 - _supply[i] = 0;
10.889 - }
10.890 - _supply[_node_id[_psource]] = _pstflow;
10.891 - _supply[_node_id[_ptarget]] = -_pstflow;
10.892 - }
10.893 - if (!valid_supply) return false;
10.894 -
10.895 - // Infinite capacity value
10.896 - Flow inf_cap =
10.897 - std::numeric_limits<Flow>::has_infinity ?
10.898 - std::numeric_limits<Flow>::infinity() :
10.899 - std::numeric_limits<Flow>::max();
10.900 -
10.901 - // Initialize artifical cost
10.902 - Cost art_cost;
10.903 - if (std::numeric_limits<Cost>::is_exact) {
10.904 - art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
10.905 - } else {
10.906 - art_cost = std::numeric_limits<Cost>::min();
10.907 - for (int i = 0; i != _arc_num; ++i) {
10.908 - if (_cost[i] > art_cost) art_cost = _cost[i];
10.909 - }
10.910 - art_cost = (art_cost + 1) * _node_num;
10.911 + ART_COST = (ART_COST + 1) * _node_num;
10.912 }
10.913
10.914 - // Run Circulation to check if a feasible solution exists
10.915 - typedef ConstMap<Arc, Flow> ConstArcMap;
10.916 - ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
10.917 - FlowNodeMap *csup = NULL;
10.918 - bool local_csup = false;
10.919 - if (_psupply) {
10.920 - csup = _psupply;
10.921 - } else {
10.922 - csup = new FlowNodeMap(_graph, 0);
10.923 - (*csup)[_psource] = _pstflow;
10.924 - (*csup)[_ptarget] = -_pstflow;
10.925 - local_csup = true;
10.926 + // Initialize arc maps
10.927 + for (int i = 0; i != _arc_num; ++i) {
10.928 + _flow[i] = 0;
10.929 + _state[i] = STATE_LOWER;
10.930 }
10.931 - bool circ_result = false;
10.932 - if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
10.933 - // GEQ problem type
10.934 - if (_plower) {
10.935 - if (_pupper) {
10.936 - Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
10.937 - circ(_graph, *_plower, *_pupper, *csup);
10.938 - circ_result = circ.run();
10.939 - } else {
10.940 - Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
10.941 - circ(_graph, *_plower, inf_arc_map, *csup);
10.942 - circ_result = circ.run();
10.943 - }
10.944 - } else {
10.945 - if (_pupper) {
10.946 - Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
10.947 - circ(_graph, zero_arc_map, *_pupper, *csup);
10.948 - circ_result = circ.run();
10.949 - } else {
10.950 - Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
10.951 - circ(_graph, zero_arc_map, inf_arc_map, *csup);
10.952 - circ_result = circ.run();
10.953 - }
10.954 - }
10.955 - } else {
10.956 - // LEQ problem type
10.957 - typedef ReverseDigraph<const GR> RevGraph;
10.958 - typedef NegMap<FlowNodeMap> NegNodeMap;
10.959 - RevGraph rgraph(_graph);
10.960 - NegNodeMap neg_csup(*csup);
10.961 - if (_plower) {
10.962 - if (_pupper) {
10.963 - Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
10.964 - circ(rgraph, *_plower, *_pupper, neg_csup);
10.965 - circ_result = circ.run();
10.966 - } else {
10.967 - Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
10.968 - circ(rgraph, *_plower, inf_arc_map, neg_csup);
10.969 - circ_result = circ.run();
10.970 - }
10.971 - } else {
10.972 - if (_pupper) {
10.973 - Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
10.974 - circ(rgraph, zero_arc_map, *_pupper, neg_csup);
10.975 - circ_result = circ.run();
10.976 - } else {
10.977 - Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
10.978 - circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
10.979 - circ_result = circ.run();
10.980 - }
10.981 - }
10.982 - }
10.983 - if (local_csup) delete csup;
10.984 - if (!circ_result) return false;
10.985 -
10.986 +
10.987 // Set data for the artificial root node
10.988 _root = _node_num;
10.989 _parent[_root] = -1;
10.990 _pred[_root] = -1;
10.991 _thread[_root] = 0;
10.992 _rev_thread[0] = _root;
10.993 - _succ_num[_root] = all_node_num;
10.994 + _succ_num[_root] = _node_num + 1;
10.995 _last_succ[_root] = _root - 1;
10.996 - _supply[_root] = -sum_supply;
10.997 - if (sum_supply < 0) {
10.998 - _pi[_root] = -art_cost;
10.999 - } else {
10.1000 - _pi[_root] = art_cost;
10.1001 - }
10.1002 -
10.1003 - // Store the arcs in a mixed order
10.1004 - int k = std::max(int(std::sqrt(double(_arc_num))), 10);
10.1005 - int i = 0;
10.1006 - for (ArcIt e(_graph); e != INVALID; ++e) {
10.1007 - _arc_ref[i] = e;
10.1008 - if ((i += k) >= _arc_num) i = (i % k) + 1;
10.1009 - }
10.1010 -
10.1011 - // Initialize arc maps
10.1012 - if (_pupper && _pcost) {
10.1013 - for (int i = 0; i != _arc_num; ++i) {
10.1014 - Arc e = _arc_ref[i];
10.1015 - _source[i] = _node_id[_graph.source(e)];
10.1016 - _target[i] = _node_id[_graph.target(e)];
10.1017 - _cap[i] = (*_pupper)[e];
10.1018 - _cost[i] = (*_pcost)[e];
10.1019 - _flow[i] = 0;
10.1020 - _state[i] = STATE_LOWER;
10.1021 - }
10.1022 - } else {
10.1023 - for (int i = 0; i != _arc_num; ++i) {
10.1024 - Arc e = _arc_ref[i];
10.1025 - _source[i] = _node_id[_graph.source(e)];
10.1026 - _target[i] = _node_id[_graph.target(e)];
10.1027 - _flow[i] = 0;
10.1028 - _state[i] = STATE_LOWER;
10.1029 - }
10.1030 - if (_pupper) {
10.1031 - for (int i = 0; i != _arc_num; ++i)
10.1032 - _cap[i] = (*_pupper)[_arc_ref[i]];
10.1033 - } else {
10.1034 - for (int i = 0; i != _arc_num; ++i)
10.1035 - _cap[i] = inf_cap;
10.1036 - }
10.1037 - if (_pcost) {
10.1038 - for (int i = 0; i != _arc_num; ++i)
10.1039 - _cost[i] = (*_pcost)[_arc_ref[i]];
10.1040 - } else {
10.1041 - for (int i = 0; i != _arc_num; ++i)
10.1042 - _cost[i] = 1;
10.1043 - }
10.1044 - }
10.1045 -
10.1046 - // Remove non-zero lower bounds
10.1047 - if (_plower) {
10.1048 - for (int i = 0; i != _arc_num; ++i) {
10.1049 - Flow c = (*_plower)[_arc_ref[i]];
10.1050 - if (c != 0) {
10.1051 - _cap[i] -= c;
10.1052 - _supply[_source[i]] -= c;
10.1053 - _supply[_target[i]] += c;
10.1054 - }
10.1055 - }
10.1056 - }
10.1057 + _supply[_root] = -_sum_supply;
10.1058 + _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
10.1059
10.1060 // Add artificial arcs and initialize the spanning tree data structure
10.1061 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
10.1062 + _parent[u] = _root;
10.1063 + _pred[u] = e;
10.1064 _thread[u] = u + 1;
10.1065 _rev_thread[u + 1] = u;
10.1066 _succ_num[u] = 1;
10.1067 _last_succ[u] = u;
10.1068 - _parent[u] = _root;
10.1069 - _pred[u] = e;
10.1070 - _cost[e] = art_cost;
10.1071 - _cap[e] = inf_cap;
10.1072 + _cost[e] = ART_COST;
10.1073 + _cap[e] = INF;
10.1074 _state[e] = STATE_TREE;
10.1075 - if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
10.1076 + if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
10.1077 _flow[e] = _supply[u];
10.1078 _forward[u] = true;
10.1079 - _pi[u] = -art_cost + _pi[_root];
10.1080 + _pi[u] = -ART_COST + _pi[_root];
10.1081 } else {
10.1082 _flow[e] = -_supply[u];
10.1083 _forward[u] = false;
10.1084 - _pi[u] = art_cost + _pi[_root];
10.1085 + _pi[u] = ART_COST + _pi[_root];
10.1086 }
10.1087 }
10.1088
10.1089 @@ -1336,13 +1148,14 @@
10.1090 }
10.1091 delta = _cap[in_arc];
10.1092 int result = 0;
10.1093 - Flow d;
10.1094 + Value d;
10.1095 int e;
10.1096
10.1097 // Search the cycle along the path form the first node to the root
10.1098 for (int u = first; u != join; u = _parent[u]) {
10.1099 e = _pred[u];
10.1100 - d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
10.1101 + d = _forward[u] ?
10.1102 + _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
10.1103 if (d < delta) {
10.1104 delta = d;
10.1105 u_out = u;
10.1106 @@ -1352,7 +1165,8 @@
10.1107 // Search the cycle along the path form the second node to the root
10.1108 for (int u = second; u != join; u = _parent[u]) {
10.1109 e = _pred[u];
10.1110 - d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
10.1111 + d = _forward[u] ?
10.1112 + (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
10.1113 if (d <= delta) {
10.1114 delta = d;
10.1115 u_out = u;
10.1116 @@ -1374,7 +1188,7 @@
10.1117 void changeFlow(bool change) {
10.1118 // Augment along the cycle
10.1119 if (delta > 0) {
10.1120 - Flow val = _state[in_arc] * delta;
10.1121 + Value val = _state[in_arc] * delta;
10.1122 _flow[in_arc] += val;
10.1123 for (int u = _source[in_arc]; u != join; u = _parent[u]) {
10.1124 _flow[_pred[u]] += _forward[u] ? -val : val;
10.1125 @@ -1526,7 +1340,7 @@
10.1126 }
10.1127
10.1128 // Execute the algorithm
10.1129 - bool start(PivotRule pivot_rule) {
10.1130 + ProblemType start(PivotRule pivot_rule) {
10.1131 // Select the pivot rule implementation
10.1132 switch (pivot_rule) {
10.1133 case FIRST_ELIGIBLE:
10.1134 @@ -1540,41 +1354,55 @@
10.1135 case ALTERING_LIST:
10.1136 return start<AlteringListPivotRule>();
10.1137 }
10.1138 - return false;
10.1139 + return INFEASIBLE; // avoid warning
10.1140 }
10.1141
10.1142 template <typename PivotRuleImpl>
10.1143 - bool start() {
10.1144 + ProblemType start() {
10.1145 PivotRuleImpl pivot(*this);
10.1146
10.1147 // Execute the Network Simplex algorithm
10.1148 while (pivot.findEnteringArc()) {
10.1149 findJoinNode();
10.1150 bool change = findLeavingArc();
10.1151 + if (delta >= INF) return UNBOUNDED;
10.1152 changeFlow(change);
10.1153 if (change) {
10.1154 updateTreeStructure();
10.1155 updatePotential();
10.1156 }
10.1157 }
10.1158 -
10.1159 - // Copy flow values to _flow_map
10.1160 - if (_plower) {
10.1161 - for (int i = 0; i != _arc_num; ++i) {
10.1162 - Arc e = _arc_ref[i];
10.1163 - _flow_map->set(e, (*_plower)[e] + _flow[i]);
10.1164 - }
10.1165 - } else {
10.1166 - for (int i = 0; i != _arc_num; ++i) {
10.1167 - _flow_map->set(_arc_ref[i], _flow[i]);
10.1168 +
10.1169 + // Check feasibility
10.1170 + if (_sum_supply < 0) {
10.1171 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
10.1172 + if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
10.1173 }
10.1174 }
10.1175 - // Copy potential values to _potential_map
10.1176 - for (NodeIt n(_graph); n != INVALID; ++n) {
10.1177 - _potential_map->set(n, _pi[_node_id[n]]);
10.1178 + else if (_sum_supply > 0) {
10.1179 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
10.1180 + if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
10.1181 + }
10.1182 + }
10.1183 + else {
10.1184 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
10.1185 + if (_flow[e] != 0) return INFEASIBLE;
10.1186 + }
10.1187 }
10.1188
10.1189 - return true;
10.1190 + // Transform the solution and the supply map to the original form
10.1191 + if (_have_lower) {
10.1192 + for (int i = 0; i != _arc_num; ++i) {
10.1193 + Value c = _lower[i];
10.1194 + if (c != 0) {
10.1195 + _flow[i] += c;
10.1196 + _supply[_source[i]] += c;
10.1197 + _supply[_target[i]] -= c;
10.1198 + }
10.1199 + }
10.1200 + }
10.1201 +
10.1202 + return OPTIMAL;
10.1203 }
10.1204
10.1205 }; //class NetworkSimplex
11.1 --- a/lemon/preflow.h Wed Apr 29 16:15:29 2009 +0100
11.2 +++ b/lemon/preflow.h Wed Apr 29 19:22:14 2009 +0100
11.3 @@ -46,13 +46,13 @@
11.4 typedef CAP CapacityMap;
11.5
11.6 /// \brief The type of the flow values.
11.7 - typedef typename CapacityMap::Value Flow;
11.8 + typedef typename CapacityMap::Value Value;
11.9
11.10 /// \brief The type of the map that stores the flow values.
11.11 ///
11.12 /// The type of the map that stores the flow values.
11.13 /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
11.14 - typedef typename Digraph::template ArcMap<Flow> FlowMap;
11.15 + typedef typename Digraph::template ArcMap<Value> FlowMap;
11.16
11.17 /// \brief Instantiates a FlowMap.
11.18 ///
11.19 @@ -84,7 +84,7 @@
11.20 /// \brief The tolerance used by the algorithm
11.21 ///
11.22 /// The tolerance used by the algorithm to handle inexact computation.
11.23 - typedef lemon::Tolerance<Flow> Tolerance;
11.24 + typedef lemon::Tolerance<Value> Tolerance;
11.25
11.26 };
11.27
11.28 @@ -125,7 +125,7 @@
11.29 ///The type of the capacity map.
11.30 typedef typename Traits::CapacityMap CapacityMap;
11.31 ///The type of the flow values.
11.32 - typedef typename Traits::Flow Flow;
11.33 + typedef typename Traits::Value Value;
11.34
11.35 ///The type of the flow map.
11.36 typedef typename Traits::FlowMap FlowMap;
11.37 @@ -151,7 +151,7 @@
11.38 Elevator* _level;
11.39 bool _local_level;
11.40
11.41 - typedef typename Digraph::template NodeMap<Flow> ExcessMap;
11.42 + typedef typename Digraph::template NodeMap<Value> ExcessMap;
11.43 ExcessMap* _excess;
11.44
11.45 Tolerance _tolerance;
11.46 @@ -470,7 +470,7 @@
11.47 }
11.48
11.49 for (NodeIt n(_graph); n != INVALID; ++n) {
11.50 - Flow excess = 0;
11.51 + Value excess = 0;
11.52 for (InArcIt e(_graph, n); e != INVALID; ++e) {
11.53 excess += (*_flow)[e];
11.54 }
11.55 @@ -519,7 +519,7 @@
11.56 _level->initFinish();
11.57
11.58 for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
11.59 - Flow rem = (*_capacity)[e] - (*_flow)[e];
11.60 + Value rem = (*_capacity)[e] - (*_flow)[e];
11.61 if (_tolerance.positive(rem)) {
11.62 Node u = _graph.target(e);
11.63 if ((*_level)[u] == _level->maxLevel()) continue;
11.64 @@ -531,7 +531,7 @@
11.65 }
11.66 }
11.67 for (InArcIt e(_graph, _source); e != INVALID; ++e) {
11.68 - Flow rem = (*_flow)[e];
11.69 + Value rem = (*_flow)[e];
11.70 if (_tolerance.positive(rem)) {
11.71 Node v = _graph.source(e);
11.72 if ((*_level)[v] == _level->maxLevel()) continue;
11.73 @@ -564,11 +564,11 @@
11.74 int num = _node_num;
11.75
11.76 while (num > 0 && n != INVALID) {
11.77 - Flow excess = (*_excess)[n];
11.78 + Value excess = (*_excess)[n];
11.79 int new_level = _level->maxLevel();
11.80
11.81 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
11.82 - Flow rem = (*_capacity)[e] - (*_flow)[e];
11.83 + Value rem = (*_capacity)[e] - (*_flow)[e];
11.84 if (!_tolerance.positive(rem)) continue;
11.85 Node v = _graph.target(e);
11.86 if ((*_level)[v] < level) {
11.87 @@ -591,7 +591,7 @@
11.88 }
11.89
11.90 for (InArcIt e(_graph, n); e != INVALID; ++e) {
11.91 - Flow rem = (*_flow)[e];
11.92 + Value rem = (*_flow)[e];
11.93 if (!_tolerance.positive(rem)) continue;
11.94 Node v = _graph.source(e);
11.95 if ((*_level)[v] < level) {
11.96 @@ -637,11 +637,11 @@
11.97
11.98 num = _node_num * 20;
11.99 while (num > 0 && n != INVALID) {
11.100 - Flow excess = (*_excess)[n];
11.101 + Value excess = (*_excess)[n];
11.102 int new_level = _level->maxLevel();
11.103
11.104 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
11.105 - Flow rem = (*_capacity)[e] - (*_flow)[e];
11.106 + Value rem = (*_capacity)[e] - (*_flow)[e];
11.107 if (!_tolerance.positive(rem)) continue;
11.108 Node v = _graph.target(e);
11.109 if ((*_level)[v] < level) {
11.110 @@ -664,7 +664,7 @@
11.111 }
11.112
11.113 for (InArcIt e(_graph, n); e != INVALID; ++e) {
11.114 - Flow rem = (*_flow)[e];
11.115 + Value rem = (*_flow)[e];
11.116 if (!_tolerance.positive(rem)) continue;
11.117 Node v = _graph.source(e);
11.118 if ((*_level)[v] < level) {
11.119 @@ -778,12 +778,12 @@
11.120
11.121 Node n;
11.122 while ((n = _level->highestActive()) != INVALID) {
11.123 - Flow excess = (*_excess)[n];
11.124 + Value excess = (*_excess)[n];
11.125 int level = _level->highestActiveLevel();
11.126 int new_level = _level->maxLevel();
11.127
11.128 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
11.129 - Flow rem = (*_capacity)[e] - (*_flow)[e];
11.130 + Value rem = (*_capacity)[e] - (*_flow)[e];
11.131 if (!_tolerance.positive(rem)) continue;
11.132 Node v = _graph.target(e);
11.133 if ((*_level)[v] < level) {
11.134 @@ -806,7 +806,7 @@
11.135 }
11.136
11.137 for (InArcIt e(_graph, n); e != INVALID; ++e) {
11.138 - Flow rem = (*_flow)[e];
11.139 + Value rem = (*_flow)[e];
11.140 if (!_tolerance.positive(rem)) continue;
11.141 Node v = _graph.source(e);
11.142 if ((*_level)[v] < level) {
11.143 @@ -897,18 +897,18 @@
11.144 ///
11.145 /// \pre Either \ref run() or \ref init() must be called before
11.146 /// using this function.
11.147 - Flow flowValue() const {
11.148 + Value flowValue() const {
11.149 return (*_excess)[_target];
11.150 }
11.151
11.152 - /// \brief Returns the flow on the given arc.
11.153 + /// \brief Returns the flow value on the given arc.
11.154 ///
11.155 - /// Returns the flow on the given arc. This method can
11.156 + /// Returns the flow value on the given arc. This method can
11.157 /// be called after the second phase of the algorithm.
11.158 ///
11.159 /// \pre Either \ref run() or \ref init() must be called before
11.160 /// using this function.
11.161 - Flow flow(const Arc& arc) const {
11.162 + Value flow(const Arc& arc) const {
11.163 return (*_flow)[arc];
11.164 }
11.165
12.1 --- a/test/lp_test.cc Wed Apr 29 16:15:29 2009 +0100
12.2 +++ b/test/lp_test.cc Wed Apr 29 19:22:14 2009 +0100
12.3 @@ -21,9 +21,7 @@
12.4 #include "test_tools.h"
12.5 #include <lemon/tolerance.h>
12.6
12.7 -#ifdef HAVE_CONFIG_H
12.8 #include <lemon/config.h>
12.9 -#endif
12.10
12.11 #ifdef LEMON_HAVE_GLPK
12.12 #include <lemon/glpk.h>
13.1 --- a/test/min_cost_flow_test.cc Wed Apr 29 16:15:29 2009 +0100
13.2 +++ b/test/min_cost_flow_test.cc Wed Apr 29 19:22:14 2009 +0100
13.3 @@ -18,6 +18,7 @@
13.4
13.5 #include <iostream>
13.6 #include <fstream>
13.7 +#include <limits>
13.8
13.9 #include <lemon/list_graph.h>
13.10 #include <lemon/lgf_reader.h>
13.11 @@ -33,57 +34,57 @@
13.12
13.13 char test_lgf[] =
13.14 "@nodes\n"
13.15 - "label sup1 sup2 sup3 sup4 sup5\n"
13.16 - " 1 20 27 0 20 30\n"
13.17 - " 2 -4 0 0 -8 -3\n"
13.18 - " 3 0 0 0 0 0\n"
13.19 - " 4 0 0 0 0 0\n"
13.20 - " 5 9 0 0 6 11\n"
13.21 - " 6 -6 0 0 -5 -6\n"
13.22 - " 7 0 0 0 0 0\n"
13.23 - " 8 0 0 0 0 3\n"
13.24 - " 9 3 0 0 0 0\n"
13.25 - " 10 -2 0 0 -7 -2\n"
13.26 - " 11 0 0 0 -10 0\n"
13.27 - " 12 -20 -27 0 -30 -20\n"
13.28 - "\n"
13.29 + "label sup1 sup2 sup3 sup4 sup5 sup6\n"
13.30 + " 1 20 27 0 30 20 30\n"
13.31 + " 2 -4 0 0 0 -8 -3\n"
13.32 + " 3 0 0 0 0 0 0\n"
13.33 + " 4 0 0 0 0 0 0\n"
13.34 + " 5 9 0 0 0 6 11\n"
13.35 + " 6 -6 0 0 0 -5 -6\n"
13.36 + " 7 0 0 0 0 0 0\n"
13.37 + " 8 0 0 0 0 0 3\n"
13.38 + " 9 3 0 0 0 0 0\n"
13.39 + " 10 -2 0 0 0 -7 -2\n"
13.40 + " 11 0 0 0 0 -10 0\n"
13.41 + " 12 -20 -27 0 -30 -30 -20\n"
13.42 + "\n"
13.43 "@arcs\n"
13.44 - " cost cap low1 low2\n"
13.45 - " 1 2 70 11 0 8\n"
13.46 - " 1 3 150 3 0 1\n"
13.47 - " 1 4 80 15 0 2\n"
13.48 - " 2 8 80 12 0 0\n"
13.49 - " 3 5 140 5 0 3\n"
13.50 - " 4 6 60 10 0 1\n"
13.51 - " 4 7 80 2 0 0\n"
13.52 - " 4 8 110 3 0 0\n"
13.53 - " 5 7 60 14 0 0\n"
13.54 - " 5 11 120 12 0 0\n"
13.55 - " 6 3 0 3 0 0\n"
13.56 - " 6 9 140 4 0 0\n"
13.57 - " 6 10 90 8 0 0\n"
13.58 - " 7 1 30 5 0 0\n"
13.59 - " 8 12 60 16 0 4\n"
13.60 - " 9 12 50 6 0 0\n"
13.61 - "10 12 70 13 0 5\n"
13.62 - "10 2 100 7 0 0\n"
13.63 - "10 7 60 10 0 0\n"
13.64 - "11 10 20 14 0 6\n"
13.65 - "12 11 30 10 0 0\n"
13.66 + " cost cap low1 low2 low3\n"
13.67 + " 1 2 70 11 0 8 8\n"
13.68 + " 1 3 150 3 0 1 0\n"
13.69 + " 1 4 80 15 0 2 2\n"
13.70 + " 2 8 80 12 0 0 0\n"
13.71 + " 3 5 140 5 0 3 1\n"
13.72 + " 4 6 60 10 0 1 0\n"
13.73 + " 4 7 80 2 0 0 0\n"
13.74 + " 4 8 110 3 0 0 0\n"
13.75 + " 5 7 60 14 0 0 0\n"
13.76 + " 5 11 120 12 0 0 0\n"
13.77 + " 6 3 0 3 0 0 0\n"
13.78 + " 6 9 140 4 0 0 0\n"
13.79 + " 6 10 90 8 0 0 0\n"
13.80 + " 7 1 30 5 0 0 -5\n"
13.81 + " 8 12 60 16 0 4 3\n"
13.82 + " 9 12 50 6 0 0 0\n"
13.83 + "10 12 70 13 0 5 2\n"
13.84 + "10 2 100 7 0 0 0\n"
13.85 + "10 7 60 10 0 0 -3\n"
13.86 + "11 10 20 14 0 6 -20\n"
13.87 + "12 11 30 10 0 0 -10\n"
13.88 "\n"
13.89 "@attributes\n"
13.90 "source 1\n"
13.91 "target 12\n";
13.92
13.93
13.94 -enum ProblemType {
13.95 +enum SupplyType {
13.96 EQ,
13.97 GEQ,
13.98 LEQ
13.99 };
13.100
13.101 // Check the interface of an MCF algorithm
13.102 -template <typename GR, typename Flow, typename Cost>
13.103 +template <typename GR, typename Value, typename Cost>
13.104 class McfClassConcept
13.105 {
13.106 public:
13.107 @@ -94,53 +95,46 @@
13.108 checkConcept<concepts::Digraph, GR>();
13.109
13.110 MCF mcf(g);
13.111 + const MCF& const_mcf = mcf;
13.112
13.113 b = mcf.reset()
13.114 .lowerMap(lower)
13.115 .upperMap(upper)
13.116 - .capacityMap(upper)
13.117 - .boundMaps(lower, upper)
13.118 .costMap(cost)
13.119 .supplyMap(sup)
13.120 .stSupply(n, n, k)
13.121 - .flowMap(flow)
13.122 - .potentialMap(pot)
13.123 .run();
13.124 -
13.125 - const MCF& const_mcf = mcf;
13.126
13.127 - const typename MCF::FlowMap &fm = const_mcf.flowMap();
13.128 - const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
13.129 -
13.130 - v = const_mcf.totalCost();
13.131 - double x = const_mcf.template totalCost<double>();
13.132 + c = const_mcf.totalCost();
13.133 + x = const_mcf.template totalCost<double>();
13.134 v = const_mcf.flow(a);
13.135 - v = const_mcf.potential(n);
13.136 -
13.137 - ignore_unused_variable_warning(fm);
13.138 - ignore_unused_variable_warning(pm);
13.139 - ignore_unused_variable_warning(x);
13.140 + c = const_mcf.potential(n);
13.141 + const_mcf.flowMap(fm);
13.142 + const_mcf.potentialMap(pm);
13.143 }
13.144
13.145 typedef typename GR::Node Node;
13.146 typedef typename GR::Arc Arc;
13.147 - typedef concepts::ReadMap<Node, Flow> NM;
13.148 - typedef concepts::ReadMap<Arc, Flow> FAM;
13.149 + typedef concepts::ReadMap<Node, Value> NM;
13.150 + typedef concepts::ReadMap<Arc, Value> VAM;
13.151 typedef concepts::ReadMap<Arc, Cost> CAM;
13.152 + typedef concepts::WriteMap<Arc, Value> FlowMap;
13.153 + typedef concepts::WriteMap<Node, Cost> PotMap;
13.154
13.155 const GR &g;
13.156 - const FAM &lower;
13.157 - const FAM &upper;
13.158 + const VAM &lower;
13.159 + const VAM &upper;
13.160 const CAM &cost;
13.161 const NM ⊃
13.162 const Node &n;
13.163 const Arc &a;
13.164 - const Flow &k;
13.165 - Flow v;
13.166 + const Value &k;
13.167 + FlowMap fm;
13.168 + PotMap pm;
13.169 bool b;
13.170 -
13.171 - typename MCF::FlowMap &flow;
13.172 - typename MCF::PotentialMap &pot;
13.173 + double x;
13.174 + typename MCF::Value v;
13.175 + typename MCF::Cost c;
13.176 };
13.177
13.178 };
13.179 @@ -151,7 +145,7 @@
13.180 typename SM, typename FM >
13.181 bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
13.182 const SM& supply, const FM& flow,
13.183 - ProblemType type = EQ )
13.184 + SupplyType type = EQ )
13.185 {
13.186 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
13.187
13.188 @@ -208,21 +202,25 @@
13.189 // Run a minimum cost flow algorithm and check the results
13.190 template < typename MCF, typename GR,
13.191 typename LM, typename UM,
13.192 - typename CM, typename SM >
13.193 -void checkMcf( const MCF& mcf, bool mcf_result,
13.194 + typename CM, typename SM,
13.195 + typename PT >
13.196 +void checkMcf( const MCF& mcf, PT mcf_result,
13.197 const GR& gr, const LM& lower, const UM& upper,
13.198 const CM& cost, const SM& supply,
13.199 - bool result, typename CM::Value total,
13.200 + PT result, bool optimal, typename CM::Value total,
13.201 const std::string &test_id = "",
13.202 - ProblemType type = EQ )
13.203 + SupplyType type = EQ )
13.204 {
13.205 check(mcf_result == result, "Wrong result " + test_id);
13.206 - if (result) {
13.207 - check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
13.208 + if (optimal) {
13.209 + typename GR::template ArcMap<typename SM::Value> flow(gr);
13.210 + typename GR::template NodeMap<typename CM::Value> pi(gr);
13.211 + mcf.flowMap(flow);
13.212 + mcf.potentialMap(pi);
13.213 + check(checkFlow(gr, lower, upper, supply, flow, type),
13.214 "The flow is not feasible " + test_id);
13.215 check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
13.216 - check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
13.217 - mcf.potentialMap()),
13.218 + check(checkPotential(gr, lower, upper, cost, supply, flow, pi),
13.219 "Wrong potentials " + test_id);
13.220 }
13.221 }
13.222 @@ -231,11 +229,13 @@
13.223 {
13.224 // Check the interfaces
13.225 {
13.226 - typedef int Flow;
13.227 - typedef int Cost;
13.228 typedef concepts::Digraph GR;
13.229 - checkConcept< McfClassConcept<GR, Flow, Cost>,
13.230 - NetworkSimplex<GR, Flow, Cost> >();
13.231 + checkConcept< McfClassConcept<GR, int, int>,
13.232 + NetworkSimplex<GR> >();
13.233 + checkConcept< McfClassConcept<GR, double, double>,
13.234 + NetworkSimplex<GR, double> >();
13.235 + checkConcept< McfClassConcept<GR, int, double>,
13.236 + NetworkSimplex<GR, int, double> >();
13.237 }
13.238
13.239 // Run various MCF tests
13.240 @@ -244,8 +244,8 @@
13.241
13.242 // Read the test digraph
13.243 Digraph gr;
13.244 - Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
13.245 - Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
13.246 + Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
13.247 + Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
13.248 ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
13.249 Node v, w;
13.250
13.251 @@ -255,14 +255,56 @@
13.252 .arcMap("cap", u)
13.253 .arcMap("low1", l1)
13.254 .arcMap("low2", l2)
13.255 + .arcMap("low3", l3)
13.256 .nodeMap("sup1", s1)
13.257 .nodeMap("sup2", s2)
13.258 .nodeMap("sup3", s3)
13.259 .nodeMap("sup4", s4)
13.260 .nodeMap("sup5", s5)
13.261 + .nodeMap("sup6", s6)
13.262 .node("source", v)
13.263 .node("target", w)
13.264 .run();
13.265 +
13.266 + // Build a test digraph for testing negative costs
13.267 + Digraph ngr;
13.268 + Node n1 = ngr.addNode();
13.269 + Node n2 = ngr.addNode();
13.270 + Node n3 = ngr.addNode();
13.271 + Node n4 = ngr.addNode();
13.272 + Node n5 = ngr.addNode();
13.273 + Node n6 = ngr.addNode();
13.274 + Node n7 = ngr.addNode();
13.275 +
13.276 + Arc a1 = ngr.addArc(n1, n2);
13.277 + Arc a2 = ngr.addArc(n1, n3);
13.278 + Arc a3 = ngr.addArc(n2, n4);
13.279 + Arc a4 = ngr.addArc(n3, n4);
13.280 + Arc a5 = ngr.addArc(n3, n2);
13.281 + Arc a6 = ngr.addArc(n5, n3);
13.282 + Arc a7 = ngr.addArc(n5, n6);
13.283 + Arc a8 = ngr.addArc(n6, n7);
13.284 + Arc a9 = ngr.addArc(n7, n5);
13.285 +
13.286 + Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
13.287 + ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
13.288 + Digraph::NodeMap<int> ns(ngr, 0);
13.289 +
13.290 + nl2[a7] = 1000;
13.291 + nl2[a8] = -1000;
13.292 +
13.293 + ns[n1] = 100;
13.294 + ns[n4] = -100;
13.295 +
13.296 + nc[a1] = 100;
13.297 + nc[a2] = 30;
13.298 + nc[a3] = 20;
13.299 + nc[a4] = 80;
13.300 + nc[a5] = 50;
13.301 + nc[a6] = 10;
13.302 + nc[a7] = 80;
13.303 + nc[a8] = 30;
13.304 + nc[a9] = -120;
13.305
13.306 // A. Test NetworkSimplex with the default pivot rule
13.307 {
13.308 @@ -271,63 +313,77 @@
13.309 // Check the equality form
13.310 mcf.upperMap(u).costMap(c);
13.311 checkMcf(mcf, mcf.supplyMap(s1).run(),
13.312 - gr, l1, u, c, s1, true, 5240, "#A1");
13.313 + gr, l1, u, c, s1, mcf.OPTIMAL, true, 5240, "#A1");
13.314 checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
13.315 - gr, l1, u, c, s2, true, 7620, "#A2");
13.316 + gr, l1, u, c, s2, mcf.OPTIMAL, true, 7620, "#A2");
13.317 mcf.lowerMap(l2);
13.318 checkMcf(mcf, mcf.supplyMap(s1).run(),
13.319 - gr, l2, u, c, s1, true, 5970, "#A3");
13.320 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#A3");
13.321 checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
13.322 - gr, l2, u, c, s2, true, 8010, "#A4");
13.323 + gr, l2, u, c, s2, mcf.OPTIMAL, true, 8010, "#A4");
13.324 mcf.reset();
13.325 checkMcf(mcf, mcf.supplyMap(s1).run(),
13.326 - gr, l1, cu, cc, s1, true, 74, "#A5");
13.327 + gr, l1, cu, cc, s1, mcf.OPTIMAL, true, 74, "#A5");
13.328 checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
13.329 - gr, l2, cu, cc, s2, true, 94, "#A6");
13.330 + gr, l2, cu, cc, s2, mcf.OPTIMAL, true, 94, "#A6");
13.331 mcf.reset();
13.332 checkMcf(mcf, mcf.run(),
13.333 - gr, l1, cu, cc, s3, true, 0, "#A7");
13.334 - checkMcf(mcf, mcf.boundMaps(l2, u).run(),
13.335 - gr, l2, u, cc, s3, false, 0, "#A8");
13.336 + gr, l1, cu, cc, s3, mcf.OPTIMAL, true, 0, "#A7");
13.337 + checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
13.338 + gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
13.339 + mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
13.340 + checkMcf(mcf, mcf.run(),
13.341 + gr, l3, u, c, s4, mcf.OPTIMAL, true, 6360, "#A9");
13.342
13.343 // Check the GEQ form
13.344 - mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
13.345 + mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
13.346 checkMcf(mcf, mcf.run(),
13.347 - gr, l1, u, c, s4, true, 3530, "#A9", GEQ);
13.348 - mcf.problemType(mcf.GEQ);
13.349 + gr, l1, u, c, s5, mcf.OPTIMAL, true, 3530, "#A10", GEQ);
13.350 + mcf.supplyType(mcf.GEQ);
13.351 checkMcf(mcf, mcf.lowerMap(l2).run(),
13.352 - gr, l2, u, c, s4, true, 4540, "#A10", GEQ);
13.353 - mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
13.354 + gr, l2, u, c, s5, mcf.OPTIMAL, true, 4540, "#A11", GEQ);
13.355 + mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
13.356 checkMcf(mcf, mcf.run(),
13.357 - gr, l2, u, c, s5, false, 0, "#A11", GEQ);
13.358 + gr, l2, u, c, s6, mcf.INFEASIBLE, false, 0, "#A12", GEQ);
13.359
13.360 // Check the LEQ form
13.361 - mcf.reset().problemType(mcf.LEQ);
13.362 - mcf.upperMap(u).costMap(c).supplyMap(s5);
13.363 + mcf.reset().supplyType(mcf.LEQ);
13.364 + mcf.upperMap(u).costMap(c).supplyMap(s6);
13.365 checkMcf(mcf, mcf.run(),
13.366 - gr, l1, u, c, s5, true, 5080, "#A12", LEQ);
13.367 + gr, l1, u, c, s6, mcf.OPTIMAL, true, 5080, "#A13", LEQ);
13.368 checkMcf(mcf, mcf.lowerMap(l2).run(),
13.369 - gr, l2, u, c, s5, true, 5930, "#A13", LEQ);
13.370 - mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
13.371 + gr, l2, u, c, s6, mcf.OPTIMAL, true, 5930, "#A14", LEQ);
13.372 + mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
13.373 checkMcf(mcf, mcf.run(),
13.374 - gr, l2, u, c, s4, false, 0, "#A14", LEQ);
13.375 + gr, l2, u, c, s5, mcf.INFEASIBLE, false, 0, "#A15", LEQ);
13.376 +
13.377 + // Check negative costs
13.378 + NetworkSimplex<Digraph> nmcf(ngr);
13.379 + nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
13.380 + checkMcf(nmcf, nmcf.run(),
13.381 + ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false, 0, "#A16");
13.382 + checkMcf(nmcf, nmcf.upperMap(nu2).run(),
13.383 + ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true, -40000, "#A17");
13.384 + nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
13.385 + checkMcf(nmcf, nmcf.run(),
13.386 + ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false, 0, "#A18");
13.387 }
13.388
13.389 // B. Test NetworkSimplex with each pivot rule
13.390 {
13.391 NetworkSimplex<Digraph> mcf(gr);
13.392 - mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
13.393 + mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
13.394
13.395 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
13.396 - gr, l2, u, c, s1, true, 5970, "#B1");
13.397 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B1");
13.398 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
13.399 - gr, l2, u, c, s1, true, 5970, "#B2");
13.400 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B2");
13.401 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
13.402 - gr, l2, u, c, s1, true, 5970, "#B3");
13.403 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B3");
13.404 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
13.405 - gr, l2, u, c, s1, true, 5970, "#B4");
13.406 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B4");
13.407 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
13.408 - gr, l2, u, c, s1, true, 5970, "#B5");
13.409 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B5");
13.410 }
13.411
13.412 return 0;
14.1 --- a/test/mip_test.cc Wed Apr 29 16:15:29 2009 +0100
14.2 +++ b/test/mip_test.cc Wed Apr 29 19:22:14 2009 +0100
14.3 @@ -18,9 +18,7 @@
14.4
14.5 #include "test_tools.h"
14.6
14.7 -#ifdef HAVE_CONFIG_H
14.8 #include <lemon/config.h>
14.9 -#endif
14.10
14.11 #ifdef LEMON_HAVE_CPLEX
14.12 #include <lemon/cplex.h>
15.1 --- a/tools/dimacs-solver.cc Wed Apr 29 16:15:29 2009 +0100
15.2 +++ b/tools/dimacs-solver.cc Wed Apr 29 19:22:14 2009 +0100
15.3 @@ -119,8 +119,8 @@
15.4
15.5 ti.restart();
15.6 NetworkSimplex<Digraph, Value> ns(g);
15.7 - ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
15.8 - if (sum_sup > 0) ns.problemType(ns.LEQ);
15.9 + ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup);
15.10 + if (sum_sup > 0) ns.supplyType(ns.LEQ);
15.11 if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';
15.12 ti.restart();
15.13 bool res = ns.run();