Changes in / [646:e01957e96c67:645:cb8270a98660] in lemon-main
- Files:
-
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
CMakeLists.txt
r631 r627 18 18 FIND_PACKAGE(COIN) 19 19 20 ADD_DEFINITIONS(-DHAVE_CONFIG_H) 21 20 22 IF(MSVC) 21 23 SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4250 /wd4355 /wd4800 /wd4996") … … 26 28 # C4996: 'function': was declared deprecated 27 29 ENDIF(MSVC) 30 31 ADD_DEFINITIONS(-DHAVE_CONFIG_H) 28 32 29 33 INCLUDE(CheckTypeSize) -
Makefile.am
r629 r567 12 12 m4/lx_check_glpk.m4 \ 13 13 m4/lx_check_soplex.m4 \ 14 m4/lx_check_coin.m4 \ 14 m4/lx_check_clp.m4 \ 15 m4/lx_check_cbc.m4 \ 15 16 CMakeLists.txt \ 16 17 cmake/FindGhostscript.cmake \ 17 cmake/FindCPLEX.cmake \18 18 cmake/FindGLPK.cmake \ 19 cmake/FindCOIN.cmake \20 19 cmake/version.cmake.in \ 21 20 cmake/version.cmake \ -
cmake/FindCOIN.cmake
r634 r627 2 2 3 3 FIND_PATH(COIN_INCLUDE_DIR coin/CoinUtilsConfig.h 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 ) 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) 46 26 47 27 INCLUDE(FindPackageHandleStandardArgs) -
cmake/FindCPLEX.cmake
r636 r627 1 SET(CPLEX_ROOT_DIR "" CACHE PATH "CPLEX root directory")2 3 1 FIND_PATH(CPLEX_INCLUDE_DIR 4 2 ilcplex/cplex.h 5 PATHS "C:/ILOG/CPLEX91/include" 6 PATHS "/opt/ilog/cplex91/include" 7 HINTS ${CPLEX_ROOT_DIR}/include 8 ) 3 PATHS "C:/ILOG/CPLEX91/include") 4 9 5 FIND_LIBRARY(CPLEX_LIBRARY 10 cplex91 11 PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda" 12 PATHS "/opt/ilog/cplex91/bin" 13 HINTS ${CPLEX_ROOT_DIR}/bin 14 ) 6 NAMES cplex91 7 PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda") 15 8 16 9 INCLUDE(FindPackageHandleStandardArgs) … … 19 12 FIND_PATH(CPLEX_BIN_DIR 20 13 cplex91.dll 21 PATHS "C:/ILOG/CPLEX91/bin/x86_win32" 22 ) 14 PATHS "C:/ILOG/CPLEX91/bin/x86_win32") 23 15 24 16 IF(CPLEX_FOUND) 25 17 SET(CPLEX_INCLUDE_DIRS ${CPLEX_INCLUDE_DIR}) 26 18 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")30 19 ENDIF(CPLEX_FOUND) 31 20 -
cmake/FindGLPK.cmake
r638 r627 1 SET(GLPK_ROOT_DIR "" CACHE PATH "GLPK root directory")2 3 1 SET(GLPK_REGKEY "[HKEY_LOCAL_MACHINE\\SOFTWARE\\GnuWin32\\Glpk;InstallPath]") 4 2 GET_FILENAME_COMPONENT(GLPK_ROOT_PATH ${GLPK_REGKEY} ABSOLUTE) … … 6 4 FIND_PATH(GLPK_INCLUDE_DIR 7 5 glpk.h 8 PATHS ${GLPK_REGKEY}/include 9 HINTS ${GLPK_ROOT_DIR}/include 10 ) 6 PATHS ${GLPK_REGKEY}/include) 7 11 8 FIND_LIBRARY(GLPK_LIBRARY 12 glpk 13 PATHS ${GLPK_REGKEY}/lib 14 HINTS ${GLPK_ROOT_DIR}/lib 15 ) 16 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) 9 NAMES glpk 10 PATHS ${GLPK_REGKEY}/lib) 45 11 46 12 INCLUDE(FindPackageHandleStandardArgs) 47 FIND_PACKAGE_HANDLE_STANDARD_ARGS(GLPK DEFAULT_MSG GLPK_LIBRARY GLPK_INCLUDE_DIR GLPK_PROPER_VERSION_FOUND)13 FIND_PACKAGE_HANDLE_STANDARD_ARGS(GLPK DEFAULT_MSG GLPK_LIBRARY GLPK_INCLUDE_DIR) 48 14 49 15 IF(GLPK_FOUND) -
doc/groups.dox
r640 r611 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, \f$lower: A\rightarrow\mathbf{Z}\f$,356 \f$ upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and355 Formally, let \f$G=(V,A)\f$ be a digraph, 356 \f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and 357 357 upper bounds for the flow values on the arcs, for which 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 flow360 on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the358 \f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$. 359 \f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow 360 on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the 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} \f$ solution365 A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution 366 366 of the following optimization problem. 367 367 … … 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} \f$ feasible solution of the problem407 An \f$f: A\rightarrow\mathbf{Z}^+_0\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 … … 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$ nodes:416 - For all \f$u\in V\f$: 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 potential function\f$\pi\f$, i.e.421 \f$uv\in A\f$ with respect to the node potentials \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 -
lemon/Makefile.am
r639 r627 16 16 lemon/bits/windows.cc 17 17 18 nodist_lemon_HEADERS = lemon/config.h 19 18 20 19 lemon_libemon_la_CXXFLAGS = \ 21 20 $(AM_CXXFLAGS) \ … … 59 58 lemon/bfs.h \ 60 59 lemon/bin_heap.h \ 61 lemon/cbc.h \62 60 lemon/circulation.h \ 63 61 lemon/clp.h \ 64 62 lemon/color.h \ 65 63 lemon/concept_check.h \ 64 lemon/config.h \ 66 65 lemon/connectivity.h \ 67 66 lemon/counter.h \ -
lemon/circulation.h
r641 r622 65 65 typedef SM SupplyMap; 66 66 67 /// \brief The type of the flow and supplyvalues.68 typedef typename SupplyMap::Value Value;67 /// \brief The type of the flow values. 68 typedef typename SupplyMap::Value Flow; 69 69 70 70 /// \brief 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< Value> FlowMap;75 typedef typename Digraph::template ArcMap<Flow> FlowMap; 76 76 77 77 /// \brief Instantiates a FlowMap. … … 105 105 /// 106 106 /// The tolerance used by the algorithm to handle inexact computation. 107 typedef lemon::Tolerance< Value> Tolerance;107 typedef lemon::Tolerance<Flow> Tolerance; 108 108 109 109 }; … … 188 188 ///The type of the digraph the algorithm runs on. 189 189 typedef typename Traits::Digraph Digraph; 190 ///The type of the flow and supplyvalues.191 typedef typename Traits:: Value Value;190 ///The type of the flow values. 191 typedef typename Traits::Flow Flow; 192 192 193 193 ///The type of the lower bound map. … … 222 222 bool _local_level; 223 223 224 typedef typename Digraph::template NodeMap< Value> ExcessMap;224 typedef typename Digraph::template NodeMap<Flow> ExcessMap; 225 225 ExcessMap* _excess; 226 226 … … 531 531 (*_excess)[_g.source(e)] -= (*_lo)[e]; 532 532 } else { 533 Valuefc = -(*_excess)[_g.target(e)];533 Flow fc = -(*_excess)[_g.target(e)]; 534 534 _flow->set(e, fc); 535 535 (*_excess)[_g.target(e)] = 0; … … 564 564 int actlevel=(*_level)[act]; 565 565 int mlevel=_node_num; 566 Valueexc=(*_excess)[act];566 Flow exc=(*_excess)[act]; 567 567 568 568 for(OutArcIt e(_g,act);e!=INVALID; ++e) { 569 569 Node v = _g.target(e); 570 Valuefc=(*_up)[e]-(*_flow)[e];570 Flow fc=(*_up)[e]-(*_flow)[e]; 571 571 if(!_tol.positive(fc)) continue; 572 572 if((*_level)[v]<actlevel) { … … 592 592 for(InArcIt e(_g,act);e!=INVALID; ++e) { 593 593 Node v = _g.source(e); 594 Valuefc=(*_flow)[e]-(*_lo)[e];594 Flow fc=(*_flow)[e]-(*_lo)[e]; 595 595 if(!_tol.positive(fc)) continue; 596 596 if((*_level)[v]<actlevel) { … … 662 662 ///@{ 663 663 664 /// \brief Returns the flow valueon the given arc.665 /// 666 /// Returns the flow valueon the given arc.664 /// \brief Returns the flow on the given arc. 665 /// 666 /// Returns the flow 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 Valueflow(const Arc& arc) const {670 Flow flow(const Arc& arc) const { 671 671 return (*_flow)[arc]; 672 672 } … … 751 751 for(NodeIt n(_g);n!=INVALID;++n) 752 752 { 753 Valuedif=-(*_supply)[n];753 Flow 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]; … … 766 766 bool checkBarrier() const 767 767 { 768 Valuedelta=0;769 Value inf_cap = std::numeric_limits<Value>::has_infinity ?770 std::numeric_limits< Value>::infinity() :771 std::numeric_limits< Value>::max();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(); 772 772 for(NodeIt n(_g);n!=INVALID;++n) 773 773 if(barrier(n)) -
lemon/core.h
r639 r617 23 23 #include <algorithm> 24 24 25 #include <lemon/co nfig.h>25 #include <lemon/core.h> 26 26 #include <lemon/bits/enable_if.h> 27 27 #include <lemon/bits/traits.h> -
lemon/network_simplex.h
r643 r618 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> 33 36 34 37 namespace lemon { … … 48 51 /// In general this class is the fastest implementation available 49 52 /// in LEMON for the minimum cost flow problem. 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. 53 /// Moreover it supports both direction of the supply/demand inequality 54 /// constraints. For more information see \ref ProblemType. 57 55 /// 58 56 /// \tparam GR The digraph type the algorithm runs on. 59 /// \tparam VThe value type used for flow amounts, capacity bounds57 /// \tparam F The value type used for flow amounts, capacity bounds 60 58 /// and supply values in the algorithm. By default it is \c int. 61 59 /// \tparam C The value type used for costs and potentials in the 62 /// algorithm. By default it is the same as \c V.60 /// algorithm. By default it is the same as \c F. 63 61 /// 64 62 /// \warning Both value types must be signed and all input data must … … 68 66 /// implementations, from which the most efficient one is used 69 67 /// by default. For more information see \ref PivotRule. 70 template <typename GR, typename V = int, typename C = V>68 template <typename GR, typename F = int, typename C = F> 71 69 class NetworkSimplex 72 70 { 73 71 public: 74 72 75 /// The type of the flow amounts, capacity bounds and supply values76 typedef V Value;77 /// The type of the arc costs73 /// The flow type of the algorithm 74 typedef F Flow; 75 /// The cost type of the algorithm 78 76 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 79 88 80 89 public: 81 90 82 /// \brief Problem type constants for the \c run() function. 83 /// 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 91 /// \brief Enum type for selecting the pivot rule. 92 /// 93 /// Enum type for selecting the pivot rule for the \ref run() 94 /// function. 95 /// 96 /// \ref NetworkSimplex provides five different pivot rule 97 /// implementations that significantly affect the running time 98 /// of the algorithm. 99 /// By default \ref BLOCK_SEARCH "Block Search" is used, which 100 /// proved to be the most efficient and the most robust on various 101 /// test inputs according to our benchmark tests. 102 /// However another pivot rule can be selected using the \ref run() 103 /// function with the proper parameter. 104 enum PivotRule { 105 106 /// The First Eligible pivot rule. 107 /// The next eligible arc is selected in a wraparound fashion 108 /// in every iteration. 109 FIRST_ELIGIBLE, 110 111 /// The Best Eligible pivot rule. 112 /// The best eligible arc is selected in every iteration. 113 BEST_ELIGIBLE, 114 115 /// The Block Search pivot rule. 116 /// A specified number of arcs are examined in every iteration 117 /// in a wraparound fashion and the best eligible arc is selected 118 /// from this block. 119 BLOCK_SEARCH, 120 121 /// The Candidate List pivot rule. 122 /// In a major iteration a candidate list is built from eligible arcs 123 /// in a wraparound fashion and in the following minor iterations 124 /// the best eligible arc is selected from this list. 125 CANDIDATE_LIST, 126 127 /// The Altering Candidate List pivot rule. 128 /// It is a modified version of the Candidate List method. 129 /// It keeps only the several best eligible arcs from the former 130 /// candidate list and extends this list in every iteration. 131 ALTERING_LIST 97 132 }; 98 133 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/demand103 /// constraints of the\ref min_cost_flow "minimum cost flow problem".104 /// 105 /// The default supplytype is \c GEQ, since this form is supported134 /// \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 106 141 /// 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()142 /// algorithm as well. 143 /// The \c LEQ problem type can be selected using the \ref problemType() 109 144 /// function. 110 145 /// 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 exact116 /// formulation of theproblem is the following.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. 117 152 /** 118 153 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] … … 130 165 CARRY_SUPPLIES = GEQ, 131 166 132 /// This option means that there are <em>"less or equal"</em>133 /// supply/demand constraints in the definition, i.e. the exact134 /// formulation of theproblem is the following.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. 135 170 /** 136 171 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] … … 148 183 SATISFY_DEMANDS = LEQ 149 184 }; 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 /// 156 /// \ref NetworkSimplex provides five different pivot rule 157 /// implementations that significantly affect the running time 158 /// of the algorithm. 159 /// By default \ref BLOCK_SEARCH "Block Search" is used, which 160 /// proved to be the most efficient and the most robust on various 161 /// test inputs according to our benchmark tests. 162 /// However another pivot rule can be selected using the \ref run() 163 /// function with the proper parameter. 164 enum PivotRule { 165 166 /// The First Eligible pivot rule. 167 /// The next eligible arc is selected in a wraparound fashion 168 /// in every iteration. 169 FIRST_ELIGIBLE, 170 171 /// The Best Eligible pivot rule. 172 /// The best eligible arc is selected in every iteration. 173 BEST_ELIGIBLE, 174 175 /// The Block Search pivot rule. 176 /// A specified number of arcs are examined in every iteration 177 /// in a wraparound fashion and the best eligible arc is selected 178 /// from this block. 179 BLOCK_SEARCH, 180 181 /// The Candidate List pivot rule. 182 /// In a major iteration a candidate list is built from eligible arcs 183 /// in a wraparound fashion and in the following minor iterations 184 /// the best eligible arc is selected from this list. 185 CANDIDATE_LIST, 186 187 /// The Altering Candidate List pivot rule. 188 /// It is a modified version of the Candidate List method. 189 /// It keeps only the several best eligible arcs from the former 190 /// candidate list and extends this list in every iteration. 191 ALTERING_LIST 192 }; 193 185 194 186 private: 195 187 196 188 TEMPLATE_DIGRAPH_TYPEDEFS(GR); 189 190 typedef typename GR::template ArcMap<Flow> FlowArcMap; 191 typedef typename GR::template ArcMap<Cost> CostArcMap; 192 typedef typename GR::template NodeMap<Flow> FlowNodeMap; 197 193 198 194 typedef std::vector<Arc> ArcVector; … … 200 196 typedef std::vector<int> IntVector; 201 197 typedef std::vector<bool> BoolVector; 202 typedef std::vector< Value> ValueVector;198 typedef std::vector<Flow> FlowVector; 203 199 typedef std::vector<Cost> CostVector; 204 200 … … 218 214 219 215 // Parameters of the problem 220 bool _have_lower; 221 SupplyType _stype; 222 Value _sum_supply; 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; 223 230 224 231 // Data structures for storing the digraph 225 232 IntNodeMap _node_id; 226 IntArcMap _arc_id;233 ArcVector _arc_ref; 227 234 IntVector _source; 228 235 IntVector _target; 229 236 230 237 // Node and arc data 231 ValueVector _lower; 232 ValueVector _upper; 233 ValueVector _cap; 238 FlowVector _cap; 234 239 CostVector _cost; 235 ValueVector _supply;236 ValueVector _flow;240 FlowVector _supply; 241 FlowVector _flow; 237 242 CostVector _pi; 238 243 … … 253 258 int first, second, right, last; 254 259 int stem, par_stem, new_stem; 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; 260 Flow delta; 265 261 266 262 private: … … 664 660 /// \param graph The digraph the algorithm runs on. 665 661 NetworkSimplex(const GR& 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()) 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) 670 668 { 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 669 LEMON_ASSERT(std::numeric_limits<Flow>::is_integer && 670 std::numeric_limits<Flow>::is_signed, 671 "The flow type of NetworkSimplex must be signed integer"); 672 LEMON_ASSERT(std::numeric_limits<Cost>::is_integer && 673 std::numeric_limits<Cost>::is_signed, 674 "The cost type of NetworkSimplex must be signed integer"); 675 } 676 677 /// Destructor. 678 ~NetworkSimplex() { 679 if (_local_flow) delete _flow_map; 680 if (_local_potential) delete _potential_map; 681 } 682 683 /// \name Parameters 684 /// The parameters of the algorithm can be specified using these 685 /// functions. 686 687 /// @{ 688 689 /// \brief Set the lower bounds on the arcs. 690 /// 691 /// This function sets the lower bounds on the arcs. 692 /// If neither this function nor \ref boundMaps() is used before 693 /// calling \ref run(), the lower bounds will be set to zero 694 /// on all arcs. 695 /// 696 /// \param map An arc map storing the lower bounds. 697 /// Its \c Value type must be convertible to the \c Flow type 698 /// of the algorithm. 699 /// 700 /// \return <tt>(*this)</tt> 701 template <typename LOWER> 702 NetworkSimplex& lowerMap(const LOWER& map) { 703 delete _plower; 704 _plower = new FlowArcMap(_graph); 705 for (ArcIt a(_graph); a != INVALID; ++a) { 706 (*_plower)[a] = map[a]; 707 } 708 return *this; 709 } 710 711 /// \brief Set the upper bounds (capacities) on the arcs. 712 /// 713 /// 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. 718 /// 719 /// \param map An arc map storing the upper bounds. 720 /// Its \c Value type must be convertible to the \c Flow type 721 /// of the algorithm. 722 /// 723 /// \return <tt>(*this)</tt> 724 template<typename UPPER> 725 NetworkSimplex& upperMap(const UPPER& map) { 726 delete _pupper; 727 _pupper = new FlowArcMap(_graph); 728 for (ArcIt a(_graph); a != INVALID; ++a) { 729 (*_pupper)[a] = map[a]; 730 } 731 return *this; 732 } 733 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 /// \brief Set the costs of the arcs. 772 /// 773 /// This function sets the costs of the arcs. 774 /// If it is not used before calling \ref run(), the costs 775 /// will be set to \c 1 on all arcs. 776 /// 777 /// \param map An arc map storing the costs. 778 /// Its \c Value type must be convertible to the \c Cost type 779 /// of the algorithm. 780 /// 781 /// \return <tt>(*this)</tt> 782 template<typename COST> 783 NetworkSimplex& costMap(const COST& map) { 784 delete _pcost; 785 _pcost = new CostArcMap(_graph); 786 for (ArcIt a(_graph); a != INVALID; ++a) { 787 (*_pcost)[a] = map[a]; 788 } 789 return *this; 790 } 791 792 /// \brief Set the supply values of the nodes. 793 /// 794 /// This function sets the supply values of the nodes. 795 /// If neither this function nor \ref stSupply() is used before 796 /// calling \ref run(), the supply of each node will be set to zero. 797 /// (It makes sense only if non-zero lower bounds are given.) 798 /// 799 /// \param map A node map storing the supply values. 800 /// Its \c Value type must be convertible to the \c Flow type 801 /// of the algorithm. 802 /// 803 /// \return <tt>(*this)</tt> 804 template<typename SUP> 805 NetworkSimplex& supplyMap(const SUP& map) { 806 delete _psupply; 807 _pstsup = false; 808 _psupply = new FlowNodeMap(_graph); 809 for (NodeIt n(_graph); n != INVALID; ++n) { 810 (*_psupply)[n] = map[n]; 811 } 812 return *this; 813 } 814 815 /// \brief Set single source and target nodes and a supply value. 816 /// 817 /// This function sets a single source node and a single target node 818 /// and the required flow value. 819 /// If neither this function nor \ref supplyMap() is used before 820 /// calling \ref run(), the supply of each node will be set to zero. 821 /// (It makes sense only if non-zero lower bounds are given.) 822 /// 823 /// \param s The source node. 824 /// \param t The target node. 825 /// \param k The required amount of flow from node \c s to node \c t 826 /// (i.e. the supply of \c s and the demand of \c t). 827 /// 828 /// \return <tt>(*this)</tt> 829 NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) { 830 delete _psupply; 831 _psupply = NULL; 832 _pstsup = true; 833 _psource = s; 834 _ptarget = t; 835 _pstflow = k; 836 return *this; 837 } 838 839 /// \brief Set the problem type. 840 /// 841 /// This function sets the problem type for the algorithm. 842 /// If it is not used before calling \ref run(), the \ref GEQ problem 843 /// type will be used. 844 /// 845 /// For more information see \ref ProblemType. 846 /// 847 /// \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 _potential_map = ↦ 885 return *this; 886 } 887 888 /// @} 889 890 /// \name Execution Control 891 /// The algorithm can be executed using \ref run(). 892 893 /// @{ 894 895 /// \brief Run the algorithm. 896 /// 897 /// This function runs the algorithm. 898 /// The paramters can be specified using functions \ref lowerMap(), 899 /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(), 900 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 901 /// \ref problemType(), \ref flowMap() and \ref potentialMap(). 902 /// For example, 903 /// \code 904 /// NetworkSimplex<ListDigraph> ns(graph); 905 /// ns.boundMaps(lower, upper).costMap(cost) 906 /// .supplyMap(sup).run(); 907 /// \endcode 908 /// 909 /// This function can be called more than once. All the parameters 910 /// that have been given are kept for the next call, unless 911 /// \ref reset() is called, thus only the modified parameters 912 /// have to be set again. See \ref reset() for examples. 913 /// 914 /// \param pivot_rule The pivot rule that will be used during the 915 /// algorithm. For more information see \ref PivotRule. 916 /// 917 /// \return \c true if a feasible flow can be found. 918 bool run(PivotRule pivot_rule = BLOCK_SEARCH) { 919 return init() && start(pivot_rule); 920 } 921 922 /// \brief Reset all the parameters that have been given before. 923 /// 924 /// This function resets all the paramaters that have been given 925 /// before using functions \ref lowerMap(), \ref upperMap(), 926 /// \ref capacityMap(), \ref boundMaps(), \ref costMap(), 927 /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 928 /// \ref flowMap() and \ref potentialMap(). 929 /// 930 /// It is useful for multiple run() calls. If this function is not 931 /// used, all the parameters given before are kept for the next 932 /// \ref run() call. 933 /// 934 /// For example, 935 /// \code 936 /// NetworkSimplex<ListDigraph> ns(graph); 937 /// 938 /// // First run 939 /// ns.lowerMap(lower).capacityMap(cap).costMap(cost) 940 /// .supplyMap(sup).run(); 941 /// 942 /// // Run again with modified cost map (reset() is not called, 943 /// // so only the cost map have to be set again) 944 /// cost[e] += 100; 945 /// ns.costMap(cost).run(); 946 /// 947 /// // Run again from scratch using reset() 948 /// // (the lower bounds will be set to zero on all arcs) 949 /// ns.reset(); 950 /// ns.capacityMap(cap).costMap(cost) 951 /// .supplyMap(sup).run(); 952 /// \endcode 953 /// 954 /// \return <tt>(*this)</tt> 955 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 973 return *this; 974 } 975 976 /// @} 977 978 /// \name Query Functions 979 /// The results of the algorithm can be obtained using these 980 /// functions.\n 981 /// The \ref run() function must be called before using them. 982 983 /// @{ 984 985 /// \brief Return the total cost of the found flow. 986 /// 987 /// This function returns the total cost of the found flow. 988 /// The complexity of the function is O(e). 989 /// 990 /// \note The return type of the function can be specified as a 991 /// template parameter. For example, 992 /// \code 993 /// ns.totalCost<double>(); 994 /// \endcode 995 /// It is useful if the total cost cannot be stored in the \c Cost 996 /// type of the algorithm, which is the default return type of the 997 /// function. 998 /// 999 /// \pre \ref run() must be called before using this function. 1000 template <typename Num> 1001 Num totalCost() const { 1002 Num c = 0; 1003 if (_pcost) { 1004 for (ArcIt e(_graph); e != INVALID; ++e) 1005 c += (*_flow_map)[e] * (*_pcost)[e]; 1006 } else { 1007 for (ArcIt e(_graph); e != INVALID; ++e) 1008 c += (*_flow_map)[e]; 1009 } 1010 return c; 1011 } 1012 1013 #ifndef DOXYGEN 1014 Cost totalCost() const { 1015 return totalCost<Cost>(); 1016 } 1017 #endif 1018 1019 /// \brief Return the flow on the given arc. 1020 /// 1021 /// This function returns the flow on the given arc. 1022 /// 1023 /// \pre \ref run() must be called before using this function. 1024 Flow flow(const Arc& a) const { 1025 return (*_flow_map)[a]; 1026 } 1027 1028 /// \brief Return a const reference to the flow map. 1029 /// 1030 /// This function returns a const reference to an arc map storing 1031 /// the found flow. 1032 /// 1033 /// \pre \ref run() must be called before using this function. 1034 const FlowMap& flowMap() const { 1035 return *_flow_map; 1036 } 1037 1038 /// \brief Return the potential (dual value) of the given node. 1039 /// 1040 /// This function returns the potential (dual value) of the 1041 /// given node. 1042 /// 1043 /// \pre \ref run() must be called before using this function. 1044 Cost potential(const Node& n) const { 1045 return (*_potential_map)[n]; 1046 } 1047 1048 /// \brief Return a const reference to the potential map 1049 /// (the dual solution). 1050 /// 1051 /// This function returns a const reference to a node map storing 1052 /// the found potentials, which form the dual solution of the 1053 /// \ref min_cost_flow "minimum cost flow" problem. 1054 /// 1055 /// \pre \ref run() must be called before using this function. 1056 const PotentialMap& potentialMap() const { 1057 return *_potential_map; 1058 } 1059 1060 /// @} 1061 1062 private: 1063 1064 // Initialize internal data structures 1065 bool init() { 1066 // Initialize result maps 1067 if (!_flow_map) { 1068 _flow_map = new FlowMap(_graph); 1069 _local_flow = true; 1070 } 1071 if (!_potential_map) { 1072 _potential_map = new PotentialMap(_graph); 1073 _local_potential = true; 1074 } 1075 1076 // Initialize vectors 678 1077 _node_num = countNodes(_graph); 679 1078 _arc_num = countArcs(_graph); 680 1079 int all_node_num = _node_num + 1; 681 1080 int all_arc_num = _arc_num + _node_num; 682 1081 if (_node_num == 0) return false; 1082 1083 _arc_ref.resize(_arc_num); 683 1084 _source.resize(all_arc_num); 684 1085 _target.resize(all_arc_num); 685 1086 686 _lower.resize(all_arc_num);687 _upper.resize(all_arc_num);688 1087 _cap.resize(all_arc_num); 689 1088 _cost.resize(all_arc_num); … … 701 1100 _state.resize(all_arc_num); 702 1101 703 // Copy the graph (store the arcs in a mixed order) 704 int i = 0; 705 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 706 _node_id[n] = i; 707 } 708 int k = std::max(int(std::sqrt(double(_arc_num))), 10); 709 i = 0; 710 for (ArcIt a(_graph); a != INVALID; ++a) { 711 _arc_id[a] = i; 712 _source[i] = _node_id[_graph.source(a)]; 713 _target[i] = _node_id[_graph.target(a)]; 714 if ((i += k) >= _arc_num) i = (i % k) + 1; 715 } 716 717 // Initialize maps 718 for (int i = 0; i != _node_num; ++i) { 719 _supply[i] = 0; 720 } 721 for (int i = 0; i != _arc_num; ++i) { 722 _lower[i] = 0; 723 _upper[i] = INF; 724 _cost[i] = 1; 725 } 726 _have_lower = false; 727 _stype = GEQ; 728 } 729 730 /// \name Parameters 731 /// The parameters of the algorithm can be specified using these 732 /// functions. 733 734 /// @{ 735 736 /// \brief Set the lower bounds on the arcs. 737 /// 738 /// This function sets the lower bounds on the arcs. 739 /// If it is not used before calling \ref run(), the lower bounds 740 /// will be set to zero on all arcs. 741 /// 742 /// \param map An arc map storing the lower bounds. 743 /// Its \c Value type must be convertible to the \c Value type 744 /// of the algorithm. 745 /// 746 /// \return <tt>(*this)</tt> 747 template <typename LowerMap> 748 NetworkSimplex& lowerMap(const LowerMap& map) { 749 _have_lower = true; 750 for (ArcIt a(_graph); a != INVALID; ++a) { 751 _lower[_arc_id[a]] = map[a]; 752 } 753 return *this; 754 } 755 756 /// \brief Set the upper bounds (capacities) on the arcs. 757 /// 758 /// This function sets the upper bounds (capacities) on the 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). 762 /// 763 /// \param map An arc map storing the upper bounds. 764 /// Its \c Value type must be convertible to the \c Value type 765 /// of the algorithm. 766 /// 767 /// \return <tt>(*this)</tt> 768 template<typename UpperMap> 769 NetworkSimplex& upperMap(const UpperMap& map) { 770 for (ArcIt a(_graph); a != INVALID; ++a) { 771 _upper[_arc_id[a]] = map[a]; 772 } 773 return *this; 774 } 775 776 /// \brief Set the costs of the arcs. 777 /// 778 /// This function sets the costs of the arcs. 779 /// If it is not used before calling \ref run(), the costs 780 /// will be set to \c 1 on all arcs. 781 /// 782 /// \param map An arc map storing the costs. 783 /// Its \c Value type must be convertible to the \c Cost type 784 /// of the algorithm. 785 /// 786 /// \return <tt>(*this)</tt> 787 template<typename CostMap> 788 NetworkSimplex& costMap(const CostMap& map) { 789 for (ArcIt a(_graph); a != INVALID; ++a) { 790 _cost[_arc_id[a]] = map[a]; 791 } 792 return *this; 793 } 794 795 /// \brief Set the supply values of the nodes. 796 /// 797 /// This function sets the supply values of the nodes. 798 /// If neither this function nor \ref stSupply() is used before 799 /// calling \ref run(), the supply of each node will be set to zero. 800 /// (It makes sense only if non-zero lower bounds are given.) 801 /// 802 /// \param map A node map storing the supply values. 803 /// Its \c Value type must be convertible to the \c Value type 804 /// of the algorithm. 805 /// 806 /// \return <tt>(*this)</tt> 807 template<typename SupplyMap> 808 NetworkSimplex& supplyMap(const SupplyMap& map) { 809 for (NodeIt n(_graph); n != INVALID; ++n) { 810 _supply[_node_id[n]] = map[n]; 811 } 812 return *this; 813 } 814 815 /// \brief Set single source and target nodes and a supply value. 816 /// 817 /// This function sets a single source node and a single target node 818 /// and the required flow value. 819 /// If neither this function nor \ref supplyMap() is used before 820 /// calling \ref run(), the supply of each node will be set to zero. 821 /// (It makes sense only if non-zero lower bounds are given.) 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 /// 827 /// \param s The source node. 828 /// \param t The target node. 829 /// \param k The required amount of flow from node \c s to node \c t 830 /// (i.e. the supply of \c s and the demand of \c t). 831 /// 832 /// \return <tt>(*this)</tt> 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; 839 return *this; 840 } 841 842 /// \brief Set the type of the supply constraints. 843 /// 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 846 /// type will be used. 847 /// 848 /// For more information see \ref SupplyType. 849 /// 850 /// \return <tt>(*this)</tt> 851 NetworkSimplex& supplyType(SupplyType supply_type) { 852 _stype = supply_type; 853 return *this; 854 } 855 856 /// @} 857 858 /// \name Execution Control 859 /// The algorithm can be executed using \ref run(). 860 861 /// @{ 862 863 /// \brief Run the algorithm. 864 /// 865 /// This function runs the algorithm. 866 /// The paramters can be specified using functions \ref lowerMap(), 867 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 868 /// \ref supplyType(). 869 /// For example, 870 /// \code 871 /// NetworkSimplex<ListDigraph> ns(graph); 872 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 873 /// .supplyMap(sup).run(); 874 /// \endcode 875 /// 876 /// This function can be called more than once. All the parameters 877 /// that have been given are kept for the next call, unless 878 /// \ref reset() is called, thus only the modified parameters 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. 882 /// 883 /// \param pivot_rule The pivot rule that will be used during the 884 /// algorithm. For more information see \ref PivotRule. 885 /// 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); 898 } 899 900 /// \brief Reset all the parameters that have been given before. 901 /// 902 /// This function resets all the paramaters that have been given 903 /// before using functions \ref lowerMap(), \ref upperMap(), 904 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(). 905 /// 906 /// It is useful for multiple run() calls. If this function is not 907 /// used, all the parameters given before are kept for the next 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. 911 /// 912 /// For example, 913 /// \code 914 /// NetworkSimplex<ListDigraph> ns(graph); 915 /// 916 /// // First run 917 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 918 /// .supplyMap(sup).run(); 919 /// 920 /// // Run again with modified cost map (reset() is not called, 921 /// // so only the cost map have to be set again) 922 /// cost[e] += 100; 923 /// ns.costMap(cost).run(); 924 /// 925 /// // Run again from scratch using reset() 926 /// // (the lower bounds will be set to zero on all arcs) 927 /// ns.reset(); 928 /// ns.upperMap(capacity).costMap(cost) 929 /// .supplyMap(sup).run(); 930 /// \endcode 931 /// 932 /// \return <tt>(*this)</tt> 933 NetworkSimplex& reset() { 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; 944 return *this; 945 } 946 947 /// @} 948 949 /// \name Query Functions 950 /// The results of the algorithm can be obtained using these 951 /// functions.\n 952 /// The \ref run() function must be called before using them. 953 954 /// @{ 955 956 /// \brief Return the total cost of the found flow. 957 /// 958 /// This function returns the total cost of the found flow. 959 /// Its complexity is O(e). 960 /// 961 /// \note The return type of the function can be specified as a 962 /// template parameter. For example, 963 /// \code 964 /// ns.totalCost<double>(); 965 /// \endcode 966 /// It is useful if the total cost cannot be stored in the \c Cost 967 /// type of the algorithm, which is the default return type of the 968 /// function. 969 /// 970 /// \pre \ref run() must be called before using this function. 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]); 977 } 978 return c; 979 } 980 981 #ifndef DOXYGEN 982 Cost totalCost() const { 983 return totalCost<Cost>(); 984 } 985 #endif 986 987 /// \brief Return the flow on the given arc. 988 /// 989 /// This function returns the flow on the given arc. 990 /// 991 /// \pre \ref run() must be called before using this function. 992 Value flow(const Arc& a) const { 993 return _flow[_arc_id[a]]; 994 } 995 996 /// \brief Return the flow map (the primal solution). 997 /// 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. 1001 /// 1002 /// \pre \ref run() must be called before using this function. 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 } 1008 } 1009 1010 /// \brief Return the potential (dual value) of the given node. 1011 /// 1012 /// This function returns the potential (dual value) of the 1013 /// given node. 1014 /// 1015 /// \pre \ref run() must be called before using this function. 1016 Cost potential(const Node& n) const { 1017 return _pi[_node_id[n]]; 1018 } 1019 1020 /// \brief Return the potential map (the dual solution). 1021 /// 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. 1026 /// 1027 /// \pre \ref run() must be called before using this function. 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 } 1033 } 1034 1035 /// @} 1036 1037 private: 1038 1039 // Initialize internal data structures 1040 bool init() { 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]; 1047 } 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) { 1102 // Initialize node related data 1103 bool valid_supply = true; 1104 Flow sum_supply = 0; 1105 if (!_pstsup && !_psupply) { 1106 _pstsup = true; 1107 _psource = _ptarget = NodeIt(_graph); 1108 _pstflow = 0; 1109 } 1110 if (_psupply) { 1111 int i = 0; 1112 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 1113 _node_id[n] = i; 1114 _supply[i] = (*_psupply)[n]; 1115 sum_supply += _supply[i]; 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(); 1053 1142 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; 1143 if (_cost[i] > art_cost) art_cost = _cost[i]; 1144 } 1145 art_cost = (art_cost + 1) * _node_num; 1146 } 1147 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; 1160 } 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(); 1057 1169 } else { 1058 _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF; 1059 } 1060 _supply[_source[i]] -= c; 1061 _supply[_target[i]] += c; 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 } 1062 1184 } 1063 1185 } else { 1064 for (int i = 0; i != _arc_num; ++i) { 1065 _cap[i] = _upper[i]; 1066 } 1067 } 1068 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]; 1077 } 1078 ART_COST = (ART_COST + 1) * _node_num; 1079 } 1080 1081 // Initialize arc maps 1082 for (int i = 0; i != _arc_num; ++i) { 1083 _flow[i] = 0; 1084 _state[i] = STATE_LOWER; 1085 } 1086 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 1087 1216 // Set data for the artificial root node 1088 1217 _root = _node_num; … … 1091 1220 _thread[_root] = 0; 1092 1221 _rev_thread[0] = _root; 1093 _succ_num[_root] = _node_num + 1;1222 _succ_num[_root] = all_node_num; 1094 1223 _last_succ[_root] = _root - 1; 1095 _supply[_root] = -_sum_supply; 1096 _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST; 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 } 1097 1285 1098 1286 // Add artificial arcs and initialize the spanning tree data structure 1099 1287 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1100 _parent[u] = _root;1101 _pred[u] = e;1102 1288 _thread[u] = u + 1; 1103 1289 _rev_thread[u + 1] = u; 1104 1290 _succ_num[u] = 1; 1105 1291 _last_succ[u] = u; 1106 _cost[e] = ART_COST; 1107 _cap[e] = INF; 1292 _parent[u] = _root; 1293 _pred[u] = e; 1294 _cost[e] = art_cost; 1295 _cap[e] = inf_cap; 1108 1296 _state[e] = STATE_TREE; 1109 if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {1297 if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) { 1110 1298 _flow[e] = _supply[u]; 1111 1299 _forward[u] = true; 1112 _pi[u] = - ART_COST+ _pi[_root];1300 _pi[u] = -art_cost + _pi[_root]; 1113 1301 } else { 1114 1302 _flow[e] = -_supply[u]; 1115 1303 _forward[u] = false; 1116 _pi[u] = ART_COST+ _pi[_root];1304 _pi[u] = art_cost + _pi[_root]; 1117 1305 } 1118 1306 } … … 1149 1337 delta = _cap[in_arc]; 1150 1338 int result = 0; 1151 Valued;1339 Flow d; 1152 1340 int e; 1153 1341 … … 1155 1343 for (int u = first; u != join; u = _parent[u]) { 1156 1344 e = _pred[u]; 1157 d = _forward[u] ? 1158 _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]); 1345 d = _forward[u] ? _flow[e] : _cap[e] - _flow[e]; 1159 1346 if (d < delta) { 1160 1347 delta = d; … … 1166 1353 for (int u = second; u != join; u = _parent[u]) { 1167 1354 e = _pred[u]; 1168 d = _forward[u] ? 1169 (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e]; 1355 d = _forward[u] ? _cap[e] - _flow[e] : _flow[e]; 1170 1356 if (d <= delta) { 1171 1357 delta = d; … … 1189 1375 // Augment along the cycle 1190 1376 if (delta > 0) { 1191 Valueval = _state[in_arc] * delta;1377 Flow val = _state[in_arc] * delta; 1192 1378 _flow[in_arc] += val; 1193 1379 for (int u = _source[in_arc]; u != join; u = _parent[u]) { … … 1341 1527 1342 1528 // Execute the algorithm 1343 ProblemTypestart(PivotRule pivot_rule) {1529 bool start(PivotRule pivot_rule) { 1344 1530 // Select the pivot rule implementation 1345 1531 switch (pivot_rule) { … … 1355 1541 return start<AlteringListPivotRule>(); 1356 1542 } 1357 return INFEASIBLE; // avoid warning1543 return false; 1358 1544 } 1359 1545 1360 1546 template <typename PivotRuleImpl> 1361 ProblemTypestart() {1547 bool start() { 1362 1548 PivotRuleImpl pivot(*this); 1363 1549 … … 1366 1552 findJoinNode(); 1367 1553 bool change = findLeavingArc(); 1368 if (delta >= INF) return UNBOUNDED;1369 1554 changeFlow(change); 1370 1555 if (change) { … … 1373 1558 } 1374 1559 } 1375 1376 // Check feasibility 1377 if (_sum_supply < 0) { 1378 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1379 if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE; 1380 } 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; 1385 } 1386 } 1387 else { 1388 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1389 if (_flow[e] != 0) return INFEASIBLE; 1390 } 1391 } 1392 1393 // Transform the solution and the supply map to the original form 1394 if (_have_lower) { 1560 1561 // Copy flow values to _flow_map 1562 if (_plower) { 1395 1563 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; 1564 Arc e = _arc_ref[i]; 1565 _flow_map->set(e, (*_plower)[e] + _flow[i]); 1566 } 1567 } else { 1568 for (int i = 0; i != _arc_num; ++i) { 1569 _flow_map->set(_arc_ref[i], _flow[i]); 1570 } 1571 } 1572 // Copy potential values to _potential_map 1573 for (NodeIt n(_graph); n != INVALID; ++n) { 1574 _potential_map->set(n, _pi[_node_id[n]]); 1575 } 1576 1577 return true; 1406 1578 } 1407 1579 -
lemon/preflow.h
r641 r611 47 47 48 48 /// \brief The type of the flow values. 49 typedef typename CapacityMap::Value Value;49 typedef typename CapacityMap::Value Flow; 50 50 51 51 /// \brief The type of the map that stores the flow values. … … 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< Value> FlowMap;55 typedef typename Digraph::template ArcMap<Flow> FlowMap; 56 56 57 57 /// \brief Instantiates a FlowMap. … … 85 85 /// 86 86 /// The tolerance used by the algorithm to handle inexact computation. 87 typedef lemon::Tolerance< Value> Tolerance;87 typedef lemon::Tolerance<Flow> Tolerance; 88 88 89 89 }; … … 126 126 typedef typename Traits::CapacityMap CapacityMap; 127 127 ///The type of the flow values. 128 typedef typename Traits:: Value Value;128 typedef typename Traits::Flow Flow; 129 129 130 130 ///The type of the flow map. … … 152 152 bool _local_level; 153 153 154 typedef typename Digraph::template NodeMap< Value> ExcessMap;154 typedef typename Digraph::template NodeMap<Flow> ExcessMap; 155 155 ExcessMap* _excess; 156 156 … … 471 471 472 472 for (NodeIt n(_graph); n != INVALID; ++n) { 473 Valueexcess = 0;473 Flow excess = 0; 474 474 for (InArcIt e(_graph, n); e != INVALID; ++e) { 475 475 excess += (*_flow)[e]; … … 520 520 521 521 for (OutArcIt e(_graph, _source); e != INVALID; ++e) { 522 Valuerem = (*_capacity)[e] - (*_flow)[e];522 Flow rem = (*_capacity)[e] - (*_flow)[e]; 523 523 if (_tolerance.positive(rem)) { 524 524 Node u = _graph.target(e); … … 532 532 } 533 533 for (InArcIt e(_graph, _source); e != INVALID; ++e) { 534 Valuerem = (*_flow)[e];534 Flow rem = (*_flow)[e]; 535 535 if (_tolerance.positive(rem)) { 536 536 Node v = _graph.source(e); … … 565 565 566 566 while (num > 0 && n != INVALID) { 567 Valueexcess = (*_excess)[n];567 Flow excess = (*_excess)[n]; 568 568 int new_level = _level->maxLevel(); 569 569 570 570 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 571 Valuerem = (*_capacity)[e] - (*_flow)[e];571 Flow rem = (*_capacity)[e] - (*_flow)[e]; 572 572 if (!_tolerance.positive(rem)) continue; 573 573 Node v = _graph.target(e); … … 592 592 593 593 for (InArcIt e(_graph, n); e != INVALID; ++e) { 594 Valuerem = (*_flow)[e];594 Flow rem = (*_flow)[e]; 595 595 if (!_tolerance.positive(rem)) continue; 596 596 Node v = _graph.source(e); … … 638 638 num = _node_num * 20; 639 639 while (num > 0 && n != INVALID) { 640 Valueexcess = (*_excess)[n];640 Flow excess = (*_excess)[n]; 641 641 int new_level = _level->maxLevel(); 642 642 643 643 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 644 Valuerem = (*_capacity)[e] - (*_flow)[e];644 Flow rem = (*_capacity)[e] - (*_flow)[e]; 645 645 if (!_tolerance.positive(rem)) continue; 646 646 Node v = _graph.target(e); … … 665 665 666 666 for (InArcIt e(_graph, n); e != INVALID; ++e) { 667 Valuerem = (*_flow)[e];667 Flow rem = (*_flow)[e]; 668 668 if (!_tolerance.positive(rem)) continue; 669 669 Node v = _graph.source(e); … … 779 779 Node n; 780 780 while ((n = _level->highestActive()) != INVALID) { 781 Valueexcess = (*_excess)[n];781 Flow 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 Valuerem = (*_capacity)[e] - (*_flow)[e];786 Flow rem = (*_capacity)[e] - (*_flow)[e]; 787 787 if (!_tolerance.positive(rem)) continue; 788 788 Node v = _graph.target(e); … … 807 807 808 808 for (InArcIt e(_graph, n); e != INVALID; ++e) { 809 Valuerem = (*_flow)[e];809 Flow rem = (*_flow)[e]; 810 810 if (!_tolerance.positive(rem)) continue; 811 811 Node v = _graph.source(e); … … 898 898 /// \pre Either \ref run() or \ref init() must be called before 899 899 /// using this function. 900 ValueflowValue() const {900 Flow flowValue() const { 901 901 return (*_excess)[_target]; 902 902 } 903 903 904 /// \brief Returns the flow valueon the given arc.905 /// 906 /// Returns the flow valueon the given arc. This method can904 /// \brief Returns the flow on the given arc. 905 /// 906 /// Returns the flow 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 Valueflow(const Arc& arc) const {911 Flow flow(const Arc& arc) const { 912 912 return (*_flow)[arc]; 913 913 } -
test/lp_test.cc
r631 r627 22 22 #include <lemon/tolerance.h> 23 23 24 #ifdef HAVE_CONFIG_H 24 25 #include <lemon/config.h> 26 #endif 25 27 26 28 #ifdef LEMON_HAVE_GLPK -
test/min_cost_flow_test.cc
r642 r615 19 19 #include <iostream> 20 20 #include <fstream> 21 #include <limits>22 21 23 22 #include <lemon/list_graph.h> … … 35 34 char test_lgf[] = 36 35 "@nodes\n" 37 "label sup1 sup2 sup3 sup4 sup5 sup6\n"38 " 1 20 27 0 3020 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 06 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 03\n"46 " 9 3 0 0 0 0 0\n"47 " 10 -2 0 0 0-7 -2\n"48 " 11 0 0 0 0-10 0\n"49 " 12 -20 -27 0 -30 - 30 -20\n"50 "\n" 36 "label sup1 sup2 sup3 sup4 sup5\n" 37 " 1 20 27 0 20 30\n" 38 " 2 -4 0 0 -8 -3\n" 39 " 3 0 0 0 0 0\n" 40 " 4 0 0 0 0 0\n" 41 " 5 9 0 0 6 11\n" 42 " 6 -6 0 0 -5 -6\n" 43 " 7 0 0 0 0 0\n" 44 " 8 0 0 0 0 3\n" 45 " 9 3 0 0 0 0\n" 46 " 10 -2 0 0 -7 -2\n" 47 " 11 0 0 0 -10 0\n" 48 " 12 -20 -27 0 -30 -20\n" 49 "\n" 51 50 "@arcs\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"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" 74 73 "\n" 75 74 "@attributes\n" … … 78 77 79 78 80 enum SupplyType {79 enum ProblemType { 81 80 EQ, 82 81 GEQ, … … 85 84 86 85 // Check the interface of an MCF algorithm 87 template <typename GR, typename Value, typename Cost>86 template <typename GR, typename Flow, typename Cost> 88 87 class McfClassConcept 89 88 { … … 96 95 97 96 MCF mcf(g); 98 const MCF& const_mcf = mcf;99 97 100 98 b = mcf.reset() 101 99 .lowerMap(lower) 102 100 .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) 106 108 .run(); 107 108 c = const_mcf.totalCost(); 109 x = const_mcf.template totalCost<double>(); 109 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 double x = const_mcf.template totalCost<double>(); 110 117 v = const_mcf.flow(a); 111 c = const_mcf.potential(n); 112 const_mcf.flowMap(fm); 113 const_mcf.potentialMap(pm); 118 v = const_mcf.potential(n); 119 120 ignore_unused_variable_warning(fm); 121 ignore_unused_variable_warning(pm); 122 ignore_unused_variable_warning(x); 114 123 } 115 124 116 125 typedef typename GR::Node Node; 117 126 typedef typename GR::Arc Arc; 118 typedef concepts::ReadMap<Node, Value> NM;119 typedef concepts::ReadMap<Arc, Value> VAM;127 typedef concepts::ReadMap<Node, Flow> NM; 128 typedef concepts::ReadMap<Arc, Flow> FAM; 120 129 typedef concepts::ReadMap<Arc, Cost> CAM; 121 typedef concepts::WriteMap<Arc, Value> FlowMap;122 typedef concepts::WriteMap<Node, Cost> PotMap;123 130 124 131 const GR &g; 125 const VAM &lower;126 const VAM &upper;132 const FAM &lower; 133 const FAM &upper; 127 134 const CAM &cost; 128 135 const NM ⊃ 129 136 const Node &n; 130 137 const Arc &a; 131 const Value &k; 132 FlowMap fm; 133 PotMap pm; 138 const Flow &k; 139 Flow v; 134 140 bool b; 135 double x; 136 typename MCF:: Value v;137 typename MCF:: Cost c;141 142 typename MCF::FlowMap &flow; 143 typename MCF::PotentialMap &pot; 138 144 }; 139 145 … … 146 152 bool checkFlow( const GR& gr, const LM& lower, const UM& upper, 147 153 const SM& supply, const FM& flow, 148 SupplyType type = EQ )154 ProblemType type = EQ ) 149 155 { 150 156 TEMPLATE_DIGRAPH_TYPEDEFS(GR); … … 203 209 template < typename MCF, typename GR, 204 210 typename LM, typename UM, 205 typename CM, typename SM, 206 typename PT > 207 void checkMcf( const MCF& mcf, PT mcf_result, 211 typename CM, typename SM > 212 void checkMcf( const MCF& mcf, bool mcf_result, 208 213 const GR& gr, const LM& lower, const UM& upper, 209 214 const CM& cost, const SM& supply, 210 PT result, bool optimal, typename CM::Value total,215 bool result, typename CM::Value total, 211 216 const std::string &test_id = "", 212 SupplyType type = EQ )217 ProblemType type = EQ ) 213 218 { 214 219 check(mcf_result == result, "Wrong result " + test_id); 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), 220 if (result) { 221 check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type), 221 222 "The flow is not feasible " + test_id); 222 223 check(mcf.totalCost() == total, "The flow is not optimal " + test_id); 223 check(checkPotential(gr, lower, upper, cost, supply, flow, pi), 224 check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(), 225 mcf.potentialMap()), 224 226 "Wrong potentials " + test_id); 225 227 } … … 230 232 // Check the interfaces 231 233 { 234 typedef int Flow; 235 typedef int Cost; 232 236 typedef concepts::Digraph GR; 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> >(); 237 checkConcept< McfClassConcept<GR, Flow, Cost>, 238 NetworkSimplex<GR, Flow, Cost> >(); 239 239 } 240 240 … … 245 245 // Read the test digraph 246 246 Digraph 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);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); 249 249 ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max()); 250 250 Node v, w; … … 256 256 .arcMap("low1", l1) 257 257 .arcMap("low2", l2) 258 .arcMap("low3", l3)259 258 .nodeMap("sup1", s1) 260 259 .nodeMap("sup2", s2) … … 262 261 .nodeMap("sup4", s4) 263 262 .nodeMap("sup5", s5) 264 .nodeMap("sup6", s6)265 263 .node("source", v) 266 264 .node("target", w) 267 265 .run(); 268 269 // Build a test digraph for testing negative costs270 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 266 309 267 // A. Test NetworkSimplex with the default pivot rule … … 314 272 mcf.upperMap(u).costMap(c); 315 273 checkMcf(mcf, mcf.supplyMap(s1).run(), 316 gr, l1, u, c, s1, mcf.OPTIMAL, true,5240, "#A1");274 gr, l1, u, c, s1, true, 5240, "#A1"); 317 275 checkMcf(mcf, mcf.stSupply(v, w, 27).run(), 318 gr, l1, u, c, s2, mcf.OPTIMAL, true,7620, "#A2");276 gr, l1, u, c, s2, true, 7620, "#A2"); 319 277 mcf.lowerMap(l2); 320 278 checkMcf(mcf, mcf.supplyMap(s1).run(), 321 gr, l2, u, c, s1, mcf.OPTIMAL, true,5970, "#A3");279 gr, l2, u, c, s1, true, 5970, "#A3"); 322 280 checkMcf(mcf, mcf.stSupply(v, w, 27).run(), 323 gr, l2, u, c, s2, mcf.OPTIMAL, true,8010, "#A4");281 gr, l2, u, c, s2, true, 8010, "#A4"); 324 282 mcf.reset(); 325 283 checkMcf(mcf, mcf.supplyMap(s1).run(), 326 gr, l1, cu, cc, s1, mcf.OPTIMAL, true,74, "#A5");284 gr, l1, cu, cc, s1, true, 74, "#A5"); 327 285 checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(), 328 gr, l2, cu, cc, s2, mcf.OPTIMAL, true,94, "#A6");286 gr, l2, cu, cc, s2, true, 94, "#A6"); 329 287 mcf.reset(); 330 288 checkMcf(mcf, mcf.run(), 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"); 289 gr, l1, cu, cc, s3, true, 0, "#A7"); 290 checkMcf(mcf, mcf.boundMaps(l2, u).run(), 291 gr, l2, u, cc, s3, false, 0, "#A8"); 337 292 338 293 // Check the GEQ form 339 mcf.reset().upperMap(u).costMap(c).supplyMap(s 5);340 checkMcf(mcf, mcf.run(), 341 gr, l1, u, c, s 5, mcf.OPTIMAL, true, 3530, "#A10", GEQ);342 mcf. supplyType(mcf.GEQ);294 mcf.reset().upperMap(u).costMap(c).supplyMap(s4); 295 checkMcf(mcf, mcf.run(), 296 gr, l1, u, c, s4, true, 3530, "#A9", GEQ); 297 mcf.problemType(mcf.GEQ); 343 298 checkMcf(mcf, mcf.lowerMap(l2).run(), 344 gr, l2, u, c, s 5, mcf.OPTIMAL, true, 4540, "#A11", GEQ);345 mcf. supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);346 checkMcf(mcf, mcf.run(), 347 gr, l2, u, c, s 6, mcf.INFEASIBLE, false, 0, "#A12", GEQ);299 gr, l2, u, c, s4, true, 4540, "#A10", GEQ); 300 mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5); 301 checkMcf(mcf, mcf.run(), 302 gr, l2, u, c, s5, false, 0, "#A11", GEQ); 348 303 349 304 // Check the LEQ form 350 mcf.reset(). supplyType(mcf.LEQ);351 mcf.upperMap(u).costMap(c).supplyMap(s 6);352 checkMcf(mcf, mcf.run(), 353 gr, l1, u, c, s 6, mcf.OPTIMAL, true, 5080, "#A13", LEQ);305 mcf.reset().problemType(mcf.LEQ); 306 mcf.upperMap(u).costMap(c).supplyMap(s5); 307 checkMcf(mcf, mcf.run(), 308 gr, l1, u, c, s5, true, 5080, "#A12", LEQ); 354 309 checkMcf(mcf, mcf.lowerMap(l2).run(), 355 gr, l2, u, c, s6, mcf.OPTIMAL, true, 5930, "#A14", LEQ); 356 mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5); 357 checkMcf(mcf, mcf.run(), 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"); 310 gr, l2, u, c, s5, true, 5930, "#A13", LEQ); 311 mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4); 312 checkMcf(mcf, mcf.run(), 313 gr, l2, u, c, s4, false, 0, "#A14", LEQ); 370 314 } 371 315 … … 373 317 { 374 318 NetworkSimplex<Digraph> mcf(gr); 375 mcf.supplyMap(s1).costMap(c). upperMap(u).lowerMap(l2);319 mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2); 376 320 377 321 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE), 378 gr, l2, u, c, s1, mcf.OPTIMAL, true,5970, "#B1");322 gr, l2, u, c, s1, true, 5970, "#B1"); 379 323 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE), 380 gr, l2, u, c, s1, mcf.OPTIMAL, true,5970, "#B2");324 gr, l2, u, c, s1, true, 5970, "#B2"); 381 325 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH), 382 gr, l2, u, c, s1, mcf.OPTIMAL, true,5970, "#B3");326 gr, l2, u, c, s1, true, 5970, "#B3"); 383 327 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST), 384 gr, l2, u, c, s1, mcf.OPTIMAL, true,5970, "#B4");328 gr, l2, u, c, s1, true, 5970, "#B4"); 385 329 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST), 386 gr, l2, u, c, s1, mcf.OPTIMAL, true,5970, "#B5");330 gr, l2, u, c, s1, true, 5970, "#B5"); 387 331 } 388 332 -
test/mip_test.cc
r631 r627 19 19 #include "test_tools.h" 20 20 21 #ifdef HAVE_CONFIG_H 21 22 #include <lemon/config.h> 23 #endif 22 24 23 25 #ifdef LEMON_HAVE_CPLEX -
tools/dimacs-solver.cc
r644 r627 120 120 ti.restart(); 121 121 NetworkSimplex<Digraph, Value> ns(g); 122 ns.lowerMap(lower). upperMap(cap).costMap(cost).supplyMap(sup);123 if (sum_sup > 0) ns. supplyType(ns.LEQ);122 ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup); 123 if (sum_sup > 0) ns.problemType(ns.LEQ); 124 124 if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n'; 125 125 ti.restart();
Note: See TracChangeset
for help on using the changeset viewer.