Changes in / [692:cb8270a98660:693:e01957e96c67] in lemon
- Files:
-
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
CMakeLists.txt
r674 r678 18 18 FIND_PACKAGE(COIN) 19 19 20 ADD_DEFINITIONS(-DHAVE_CONFIG_H)21 22 20 IF(MSVC) 23 21 SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4250 /wd4355 /wd4800 /wd4996") … … 28 26 # C4996: 'function': was declared deprecated 29 27 ENDIF(MSVC) 30 31 ADD_DEFINITIONS(-DHAVE_CONFIG_H)32 28 33 29 INCLUDE(CheckTypeSize) -
Makefile.am
r614 r676 12 12 m4/lx_check_glpk.m4 \ 13 13 m4/lx_check_soplex.m4 \ 14 m4/lx_check_clp.m4 \ 15 m4/lx_check_cbc.m4 \ 14 m4/lx_check_coin.m4 \ 16 15 CMakeLists.txt \ 17 16 cmake/FindGhostscript.cmake \ 17 cmake/FindCPLEX.cmake \ 18 18 cmake/FindGLPK.cmake \ 19 cmake/FindCOIN.cmake \ 19 20 cmake/version.cmake.in \ 20 21 cmake/version.cmake \ -
cmake/FindCOIN.cmake
r674 r681 2 2 3 3 FIND_PATH(COIN_INCLUDE_DIR coin/CoinUtilsConfig.h 4 PATHS ${COIN_ROOT_DIR}/include) 5 6 FIND_LIBRARY(COIN_CBC_LIBRARY libCbc 7 PATHS ${COIN_ROOT_DIR}/lib) 8 FIND_LIBRARY(COIN_CBC_SOLVER_LIBRARY libCbcSolver 9 PATHS ${COIN_ROOT_DIR}/lib) 10 FIND_LIBRARY(COIN_CGL_LIBRARY libCgl 11 PATHS ${COIN_ROOT_DIR}/lib) 12 FIND_LIBRARY(COIN_CLP_LIBRARY libClp 13 PATHS ${COIN_ROOT_DIR}/lib) 14 FIND_LIBRARY(COIN_COIN_UTILS_LIBRARY libCoinUtils 15 PATHS ${COIN_ROOT_DIR}/lib) 16 FIND_LIBRARY(COIN_OSI_LIBRARY libOsi 17 PATHS ${COIN_ROOT_DIR}/lib) 18 FIND_LIBRARY(COIN_OSI_CBC_LIBRARY libOsiCbc 19 PATHS ${COIN_ROOT_DIR}/lib) 20 FIND_LIBRARY(COIN_OSI_CLP_LIBRARY libOsiClp 21 PATHS ${COIN_ROOT_DIR}/lib) 22 FIND_LIBRARY(COIN_OSI_VOL_LIBRARY libOsiVol 23 PATHS ${COIN_ROOT_DIR}/lib) 24 FIND_LIBRARY(COIN_VOL_LIBRARY libVol 25 PATHS ${COIN_ROOT_DIR}/lib) 4 HINTS ${COIN_ROOT_DIR}/include 5 ) 6 FIND_LIBRARY(COIN_CBC_LIBRARY 7 NAMES Cbc libCbc 8 HINTS ${COIN_ROOT_DIR}/lib 9 ) 10 FIND_LIBRARY(COIN_CBC_SOLVER_LIBRARY 11 NAMES CbcSolver libCbcSolver 12 HINTS ${COIN_ROOT_DIR}/lib 13 ) 14 FIND_LIBRARY(COIN_CGL_LIBRARY 15 NAMES Cgl libCgl 16 HINTS ${COIN_ROOT_DIR}/lib 17 ) 18 FIND_LIBRARY(COIN_CLP_LIBRARY 19 NAMES Clp libClp 20 HINTS ${COIN_ROOT_DIR}/lib 21 ) 22 FIND_LIBRARY(COIN_COIN_UTILS_LIBRARY 23 NAMES CoinUtils libCoinUtils 24 HINTS ${COIN_ROOT_DIR}/lib 25 ) 26 FIND_LIBRARY(COIN_OSI_LIBRARY 27 NAMES Osi libOsi 28 HINTS ${COIN_ROOT_DIR}/lib 29 ) 30 FIND_LIBRARY(COIN_OSI_CBC_LIBRARY 31 NAMES OsiCbc libOsiCbc 32 HINTS ${COIN_ROOT_DIR}/lib 33 ) 34 FIND_LIBRARY(COIN_OSI_CLP_LIBRARY 35 NAMES OsiClp libOsiClp 36 HINTS ${COIN_ROOT_DIR}/lib 37 ) 38 FIND_LIBRARY(COIN_OSI_VOL_LIBRARY 39 NAMES OsiVol libOsiVol 40 HINTS ${COIN_ROOT_DIR}/lib 41 ) 42 FIND_LIBRARY(COIN_VOL_LIBRARY 43 NAMES Vol libVol 44 HINTS ${COIN_ROOT_DIR}/lib 45 ) 26 46 27 47 INCLUDE(FindPackageHandleStandardArgs) -
cmake/FindCPLEX.cmake
r674 r683 1 SET(CPLEX_ROOT_DIR "" CACHE PATH "CPLEX root directory") 2 1 3 FIND_PATH(CPLEX_INCLUDE_DIR 2 4 ilcplex/cplex.h 3 PATHS "C:/ILOG/CPLEX91/include") 4 5 PATHS "C:/ILOG/CPLEX91/include" 6 PATHS "/opt/ilog/cplex91/include" 7 HINTS ${CPLEX_ROOT_DIR}/include 8 ) 5 9 FIND_LIBRARY(CPLEX_LIBRARY 6 NAMES cplex91 7 PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda") 10 cplex91 11 PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda" 12 PATHS "/opt/ilog/cplex91/bin" 13 HINTS ${CPLEX_ROOT_DIR}/bin 14 ) 8 15 9 16 INCLUDE(FindPackageHandleStandardArgs) … … 12 19 FIND_PATH(CPLEX_BIN_DIR 13 20 cplex91.dll 14 PATHS "C:/ILOG/CPLEX91/bin/x86_win32") 21 PATHS "C:/ILOG/CPLEX91/bin/x86_win32" 22 ) 15 23 16 24 IF(CPLEX_FOUND) 17 25 SET(CPLEX_INCLUDE_DIRS ${CPLEX_INCLUDE_DIR}) 18 26 SET(CPLEX_LIBRARIES ${CPLEX_LIBRARY}) 27 IF(CMAKE_SYSTEM_NAME STREQUAL "Linux") 28 SET(CPLEX_LIBRARIES "${CPLEX_LIBRARIES};m;pthread") 29 ENDIF(CMAKE_SYSTEM_NAME STREQUAL "Linux") 19 30 ENDIF(CPLEX_FOUND) 20 31 -
cmake/FindGLPK.cmake
r674 r685 1 SET(GLPK_ROOT_DIR "" CACHE PATH "GLPK root directory") 2 1 3 SET(GLPK_REGKEY "[HKEY_LOCAL_MACHINE\\SOFTWARE\\GnuWin32\\Glpk;InstallPath]") 2 4 GET_FILENAME_COMPONENT(GLPK_ROOT_PATH ${GLPK_REGKEY} ABSOLUTE) … … 4 6 FIND_PATH(GLPK_INCLUDE_DIR 5 7 glpk.h 6 PATHS ${GLPK_REGKEY}/include) 8 PATHS ${GLPK_REGKEY}/include 9 HINTS ${GLPK_ROOT_DIR}/include 10 ) 11 FIND_LIBRARY(GLPK_LIBRARY 12 glpk 13 PATHS ${GLPK_REGKEY}/lib 14 HINTS ${GLPK_ROOT_DIR}/lib 15 ) 7 16 8 FIND_LIBRARY(GLPK_LIBRARY 9 NAMES glpk 10 PATHS ${GLPK_REGKEY}/lib) 17 IF(GLPK_INCLUDE_DIR AND GLPK_LIBRARY) 18 FILE(READ ${GLPK_INCLUDE_DIR}/glpk.h GLPK_GLPK_H) 19 20 STRING(REGEX MATCH "define[ ]+GLP_MAJOR_VERSION[ ]+[0-9]+" GLPK_MAJOR_VERSION_LINE "${GLPK_GLPK_H}") 21 STRING(REGEX REPLACE "define[ ]+GLP_MAJOR_VERSION[ ]+([0-9]+)" "\\1" GLPK_VERSION_MAJOR "${GLPK_MAJOR_VERSION_LINE}") 22 23 STRING(REGEX MATCH "define[ ]+GLP_MINOR_VERSION[ ]+[0-9]+" GLPK_MINOR_VERSION_LINE "${GLPK_GLPK_H}") 24 STRING(REGEX REPLACE "define[ ]+GLP_MINOR_VERSION[ ]+([0-9]+)" "\\1" GLPK_VERSION_MINOR "${GLPK_MINOR_VERSION_LINE}") 25 26 SET(GLPK_VERSION_STRING "${GLPK_VERSION_MAJOR}.${GLPK_VERSION_MINOR}") 27 28 IF(GLPK_FIND_VERSION) 29 IF(GLPK_FIND_VERSION_COUNT GREATER 2) 30 MESSAGE(SEND_ERROR "unexpected version string") 31 ENDIF(GLPK_FIND_VERSION_COUNT GREATER 2) 32 33 MATH(EXPR GLPK_REQUESTED_VERSION "${GLPK_FIND_VERSION_MAJOR}*100 + ${GLPK_FIND_VERSION_MINOR}") 34 MATH(EXPR GLPK_FOUND_VERSION "${GLPK_VERSION_MAJOR}*100 + ${GLPK_VERSION_MINOR}") 35 36 IF(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION) 37 SET(GLPK_PROPER_VERSION_FOUND FALSE) 38 ELSE(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION) 39 SET(GLPK_PROPER_VERSION_FOUND TRUE) 40 ENDIF(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION) 41 ELSE(GLPK_FIND_VERSION) 42 SET(GLPK_PROPER_VERSION_FOUND TRUE) 43 ENDIF(GLPK_FIND_VERSION) 44 ENDIF(GLPK_INCLUDE_DIR AND GLPK_LIBRARY) 11 45 12 46 INCLUDE(FindPackageHandleStandardArgs) 13 FIND_PACKAGE_HANDLE_STANDARD_ARGS(GLPK DEFAULT_MSG GLPK_LIBRARY GLPK_INCLUDE_DIR )47 FIND_PACKAGE_HANDLE_STANDARD_ARGS(GLPK DEFAULT_MSG GLPK_LIBRARY GLPK_INCLUDE_DIR GLPK_PROPER_VERSION_FOUND) 14 48 15 49 IF(GLPK_FOUND) -
doc/groups.dox
r658 r687 353 353 in a network with capacity constraints (lower and upper bounds) 354 354 and arc costs. 355 Formally, let \f$G=(V,A)\f$ be a digraph, 356 \f$ lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and355 Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$, 356 \f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and 357 357 upper bounds for the flow values on the arcs, for which 358 \f$ 0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.359 \f$cost: A\rightarrow\mathbf{Z} ^+_0\f$ denotes the cost per unit flow360 on the arcs ,and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the358 \f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$, 359 \f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow 360 on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the 361 361 signed supply values of the nodes. 362 362 If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$ 363 363 supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with 364 364 \f$-sup(u)\f$ demand. 365 A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z} ^+_0\f$ solution365 A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\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} ^+_0\f$ feasible solution of the problem407 An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem 408 408 is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$ 409 409 node potentials the following \e complementary \e slackness optimality … … 414 414 - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$; 415 415 - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$. 416 - For all \f$u\in V\f$ :416 - For all \f$u\in V\f$ nodes: 417 417 - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$, 418 418 then \f$\pi(u)=0\f$. 419 419 420 420 Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc 421 \f$uv\in A\f$ with respect to the node potentials\f$\pi\f$, i.e.421 \f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e. 422 422 \f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f] 423 423 424 All algorithms provide dual solution (node potentials) as well 424 All algorithms provide dual solution (node potentials) as well, 425 425 if an optimal flow is found. 426 426 -
lemon/Makefile.am
r674 r686 16 16 lemon/bits/windows.cc 17 17 18 18 nodist_lemon_HEADERS = lemon/config.h 19 19 20 lemon_libemon_la_CXXFLAGS = \ 20 21 $(AM_CXXFLAGS) \ … … 58 59 lemon/bfs.h \ 59 60 lemon/bin_heap.h \ 61 lemon/cbc.h \ 60 62 lemon/circulation.h \ 61 63 lemon/clp.h \ 62 64 lemon/color.h \ 63 65 lemon/concept_check.h \ 64 lemon/config.h \65 66 lemon/connectivity.h \ 66 67 lemon/counter.h \ -
lemon/circulation.h
r669 r688 65 65 typedef SM SupplyMap; 66 66 67 /// \brief The type of the flow values.68 typedef typename SupplyMap::Value Flow;67 /// \brief The type of the flow and supply values. 68 typedef typename SupplyMap::Value Value; 69 69 70 70 /// \brief The type of the map that stores the flow values. … … 73 73 /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap" 74 74 /// concept. 75 typedef typename Digraph::template ArcMap< Flow> FlowMap;75 typedef typename Digraph::template ArcMap<Value> FlowMap; 76 76 77 77 /// \brief Instantiates a FlowMap. … … 105 105 /// 106 106 /// The tolerance used by the algorithm to handle inexact computation. 107 typedef lemon::Tolerance< Flow> Tolerance;107 typedef lemon::Tolerance<Value> Tolerance; 108 108 109 109 }; … … 188 188 ///The type of the digraph the algorithm runs on. 189 189 typedef typename Traits::Digraph Digraph; 190 ///The type of the flow values.191 typedef typename Traits:: Flow Flow;190 ///The type of the flow and supply values. 191 typedef typename Traits::Value Value; 192 192 193 193 ///The type of the lower bound map. … … 222 222 bool _local_level; 223 223 224 typedef typename Digraph::template NodeMap< Flow> ExcessMap;224 typedef typename Digraph::template NodeMap<Value> ExcessMap; 225 225 ExcessMap* _excess; 226 226 … … 531 531 (*_excess)[_g.source(e)] -= (*_lo)[e]; 532 532 } else { 533 Flowfc = -(*_excess)[_g.target(e)];533 Value 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 Flowexc=(*_excess)[act];566 Value exc=(*_excess)[act]; 567 567 568 568 for(OutArcIt e(_g,act);e!=INVALID; ++e) { 569 569 Node v = _g.target(e); 570 Flowfc=(*_up)[e]-(*_flow)[e];570 Value 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 Flowfc=(*_flow)[e]-(*_lo)[e];594 Value 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 on the given arc.665 /// 666 /// Returns the flow on the given arc.664 /// \brief Returns the flow value on the given arc. 665 /// 666 /// Returns the flow value on the given arc. 667 667 /// 668 668 /// \pre Either \ref run() or \ref init() must be called before 669 669 /// using this function. 670 Flowflow(const Arc& arc) const {670 Value flow(const Arc& arc) const { 671 671 return (*_flow)[arc]; 672 672 } … … 751 751 for(NodeIt n(_g);n!=INVALID;++n) 752 752 { 753 Flowdif=-(*_supply)[n];753 Value dif=-(*_supply)[n]; 754 754 for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e]; 755 755 for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e]; … … 766 766 bool checkBarrier() const 767 767 { 768 Flowdelta=0;769 Flow inf_cap = std::numeric_limits<Flow>::has_infinity ?770 std::numeric_limits< Flow>::infinity() :771 std::numeric_limits< Flow>::max();768 Value delta=0; 769 Value inf_cap = std::numeric_limits<Value>::has_infinity ? 770 std::numeric_limits<Value>::infinity() : 771 std::numeric_limits<Value>::max(); 772 772 for(NodeIt n(_g);n!=INVALID;++n) 773 773 if(barrier(n)) -
lemon/core.h
r664 r686 23 23 #include <algorithm> 24 24 25 #include <lemon/co re.h>25 #include <lemon/config.h> 26 26 #include <lemon/bits/enable_if.h> 27 27 #include <lemon/bits/traits.h> -
lemon/network_simplex.h
r665 r690 31 31 #include <lemon/core.h> 32 32 #include <lemon/math.h> 33 #include <lemon/maps.h>34 #include <lemon/circulation.h>35 #include <lemon/adaptors.h>36 33 37 34 namespace lemon { … … 51 48 /// In general this class is the fastest implementation available 52 49 /// in LEMON for the minimum cost flow problem. 53 /// Moreover it supports both direction of the supply/demand inequality 54 /// constraints. For more information see \ref ProblemType. 50 /// Moreover it supports both directions of the supply/demand inequality 51 /// constraints. For more information see \ref SupplyType. 52 /// 53 /// Most of the parameters of the problem (except for the digraph) 54 /// can be given using separate functions, and the algorithm can be 55 /// executed using the \ref run() function. If some parameters are not 56 /// specified, then default values will be used. 55 57 /// 56 58 /// \tparam GR The digraph type the algorithm runs on. 57 /// \tparam FThe value type used for flow amounts, capacity bounds59 /// \tparam V The value type used for flow amounts, capacity bounds 58 60 /// and supply values in the algorithm. By default it is \c int. 59 61 /// \tparam C The value type used for costs and potentials in the 60 /// algorithm. By default it is the same as \c F.62 /// algorithm. By default it is the same as \c V. 61 63 /// 62 64 /// \warning Both value types must be signed and all input data must … … 66 68 /// implementations, from which the most efficient one is used 67 69 /// by default. For more information see \ref PivotRule. 68 template <typename GR, typename F = int, typename C = F>70 template <typename GR, typename V = int, typename C = V> 69 71 class NetworkSimplex 70 72 { 71 73 public: 72 74 73 /// The flow type of the algorithm74 typedef F Flow;75 /// The cost type of the algorithm75 /// The type of the flow amounts, capacity bounds and supply values 76 typedef V Value; 77 /// The type of the arc costs 76 78 typedef C Cost; 77 #ifdef DOXYGEN78 /// The type of the flow map79 typedef GR::ArcMap<Flow> FlowMap;80 /// The type of the potential map81 typedef GR::NodeMap<Cost> PotentialMap;82 #else83 /// The type of the flow map84 typedef typename GR::template ArcMap<Flow> FlowMap;85 /// The type of the potential map86 typedef typename GR::template NodeMap<Cost> PotentialMap;87 #endif88 79 89 80 public: 90 81 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 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 132 97 }; 133 98 134 /// \brief Enum type for selecting the problem type.135 /// 136 /// Enum type for selecting the problem type, i.e. the direction of137 /// the inequalities in the supply/demand constraints of the138 /// \ref min_cost_flow "minimum cost flow problem".139 /// 140 /// The default problemtype is \c GEQ, since this form is supported99 /// \brief Constants for selecting the type of the supply constraints. 100 /// 101 /// Enum type containing constants for selecting the supply type, 102 /// i.e. the direction of the inequalities in the supply/demand 103 /// constraints of the \ref min_cost_flow "minimum cost flow problem". 104 /// 105 /// The default supply type is \c GEQ, since this form is supported 141 106 /// by other minimum cost flow algorithms and the \ref Circulation 142 /// algorithm as well.143 /// The \c LEQ problem type can be selected using the \ref problemType()107 /// algorithm, as well. 108 /// The \c LEQ problem type can be selected using the \ref supplyType() 144 109 /// function. 145 110 /// 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 the151 /// problem is the following.111 /// Note that the equality form is a special case of both supply types. 112 enum SupplyType { 113 114 /// This option means that there are <em>"greater or equal"</em> 115 /// supply/demand constraints in the definition, i.e. the exact 116 /// formulation of the problem is the following. 152 117 /** 153 118 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] … … 165 130 CARRY_SUPPLIES = GEQ, 166 131 167 /// This option means that there are "<em>less or equal</em>"168 /// constraints in the defintion, i.e. the exact formulation of the169 /// problem is the following.132 /// This option means that there are <em>"less or equal"</em> 133 /// supply/demand constraints in the definition, i.e. the exact 134 /// formulation of the problem is the following. 170 135 /** 171 136 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] … … 183 148 SATISFY_DEMANDS = LEQ 184 149 }; 185 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 186 194 private: 187 195 188 196 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;193 197 194 198 typedef std::vector<Arc> ArcVector; … … 196 200 typedef std::vector<int> IntVector; 197 201 typedef std::vector<bool> BoolVector; 198 typedef std::vector< Flow> FlowVector;202 typedef std::vector<Value> ValueVector; 199 203 typedef std::vector<Cost> CostVector; 200 204 … … 214 218 215 219 // Parameters of the problem 216 FlowArcMap *_plower; 217 FlowArcMap *_pupper; 218 CostArcMap *_pcost; 219 FlowNodeMap *_psupply; 220 bool _pstsup; 221 Node _psource, _ptarget; 222 Flow _pstflow; 223 ProblemType _ptype; 224 225 // Result maps 226 FlowMap *_flow_map; 227 PotentialMap *_potential_map; 228 bool _local_flow; 229 bool _local_potential; 220 bool _have_lower; 221 SupplyType _stype; 222 Value _sum_supply; 230 223 231 224 // Data structures for storing the digraph 232 225 IntNodeMap _node_id; 233 ArcVector _arc_ref;226 IntArcMap _arc_id; 234 227 IntVector _source; 235 228 IntVector _target; 236 229 237 230 // Node and arc data 238 FlowVector _cap; 231 ValueVector _lower; 232 ValueVector _upper; 233 ValueVector _cap; 239 234 CostVector _cost; 240 FlowVector _supply;241 FlowVector _flow;235 ValueVector _supply; 236 ValueVector _flow; 242 237 CostVector _pi; 243 238 … … 258 253 int first, second, right, last; 259 254 int stem, par_stem, new_stem; 260 Flow delta; 255 Value delta; 256 257 public: 258 259 /// \brief Constant for infinite upper bounds (capacities). 260 /// 261 /// Constant for infinite upper bounds (capacities). 262 /// It is \c std::numeric_limits<Value>::infinity() if available, 263 /// \c std::numeric_limits<Value>::max() otherwise. 264 const Value INF; 261 265 262 266 private: … … 660 664 /// \param graph The digraph the algorithm runs on. 661 665 NetworkSimplex(const GR& graph) : 662 _graph(graph), 663 _plower(NULL), _pupper(NULL), _pcost(NULL), 664 _psupply(NULL), _pstsup(false), _ptype(GEQ), 665 _flow_map(NULL), _potential_map(NULL), 666 _local_flow(false), _local_potential(false), 667 _node_id(graph) 666 _graph(graph), _node_id(graph), _arc_id(graph), 667 INF(std::numeric_limits<Value>::has_infinity ? 668 std::numeric_limits<Value>::infinity() : 669 std::numeric_limits<Value>::max()) 668 670 { 669 LEMON_ASSERT(std::numeric_limits<Flow>::is_integer && 670 std::numeric_limits<Flow>::is_signed, 671 "The flow type of NetworkSimplex must be signed integer"); 672 LEMON_ASSERT(std::numeric_limits<Cost>::is_integer && 673 std::numeric_limits<Cost>::is_signed, 674 "The cost type of NetworkSimplex must be signed integer"); 675 } 676 677 /// Destructor. 678 ~NetworkSimplex() { 679 if (_local_flow) delete _flow_map; 680 if (_local_potential) delete _potential_map; 671 // Check the value types 672 LEMON_ASSERT(std::numeric_limits<Value>::is_signed, 673 "The flow type of NetworkSimplex must be signed"); 674 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, 675 "The cost type of NetworkSimplex must be signed"); 676 677 // Resize vectors 678 _node_num = countNodes(_graph); 679 _arc_num = countArcs(_graph); 680 int all_node_num = _node_num + 1; 681 int all_arc_num = _arc_num + _node_num; 682 683 _source.resize(all_arc_num); 684 _target.resize(all_arc_num); 685 686 _lower.resize(all_arc_num); 687 _upper.resize(all_arc_num); 688 _cap.resize(all_arc_num); 689 _cost.resize(all_arc_num); 690 _supply.resize(all_node_num); 691 _flow.resize(all_arc_num); 692 _pi.resize(all_node_num); 693 694 _parent.resize(all_node_num); 695 _pred.resize(all_node_num); 696 _forward.resize(all_node_num); 697 _thread.resize(all_node_num); 698 _rev_thread.resize(all_node_num); 699 _succ_num.resize(all_node_num); 700 _last_succ.resize(all_node_num); 701 _state.resize(all_arc_num); 702 703 // Copy the graph (store the arcs in a mixed order) 704 int i = 0; 705 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 706 _node_id[n] = i; 707 } 708 int k = std::max(int(std::sqrt(double(_arc_num))), 10); 709 i = 0; 710 for (ArcIt a(_graph); a != INVALID; ++a) { 711 _arc_id[a] = i; 712 _source[i] = _node_id[_graph.source(a)]; 713 _target[i] = _node_id[_graph.target(a)]; 714 if ((i += k) >= _arc_num) i = (i % k) + 1; 715 } 716 717 // Initialize maps 718 for (int i = 0; i != _node_num; ++i) { 719 _supply[i] = 0; 720 } 721 for (int i = 0; i != _arc_num; ++i) { 722 _lower[i] = 0; 723 _upper[i] = INF; 724 _cost[i] = 1; 725 } 726 _have_lower = false; 727 _stype = GEQ; 681 728 } 682 729 … … 690 737 /// 691 738 /// This function sets the lower bounds on the arcs. 692 /// If neither this function nor \ref boundMaps() is used before 693 /// calling \ref run(), the lower bounds will be set to zero 694 /// on all arcs. 739 /// If it is not used before calling \ref run(), the lower bounds 740 /// will be set to zero on all arcs. 695 741 /// 696 742 /// \param map An arc map storing the lower bounds. 697 /// Its \c Value type must be convertible to the \c Flowtype743 /// Its \c Value type must be convertible to the \c Value type 698 744 /// of the algorithm. 699 745 /// 700 746 /// \return <tt>(*this)</tt> 701 template <typename LOWER> 702 NetworkSimplex& lowerMap(const LOWER& map) { 703 delete _plower; 704 _plower = new FlowArcMap(_graph); 747 template <typename LowerMap> 748 NetworkSimplex& lowerMap(const LowerMap& map) { 749 _have_lower = true; 705 750 for (ArcIt a(_graph); a != INVALID; ++a) { 706 (*_plower)[a] = map[a];751 _lower[_arc_id[a]] = map[a]; 707 752 } 708 753 return *this; … … 712 757 /// 713 758 /// This function sets the upper bounds (capacities) on the arcs. 714 /// If none of the functions \ref upperMap(), \ref capacityMap() 715 /// and \ref boundMaps() is used before calling \ref run(), 716 /// the upper bounds (capacities) will be set to 717 /// \c std::numeric_limits<Flow>::max() on all arcs. 759 /// If it is not used before calling \ref run(), the upper bounds 760 /// will be set to \ref INF on all arcs (i.e. the flow value will be 761 /// unbounded from above on each arc). 718 762 /// 719 763 /// \param map An arc map storing the upper bounds. 720 /// Its \c Value type must be convertible to the \c Flowtype764 /// Its \c Value type must be convertible to the \c Value type 721 765 /// of the algorithm. 722 766 /// 723 767 /// \return <tt>(*this)</tt> 724 template<typename UPPER> 725 NetworkSimplex& upperMap(const UPPER& map) { 726 delete _pupper; 727 _pupper = new FlowArcMap(_graph); 768 template<typename UpperMap> 769 NetworkSimplex& upperMap(const UpperMap& map) { 728 770 for (ArcIt a(_graph); a != INVALID; ++a) { 729 (*_pupper)[a] = map[a];771 _upper[_arc_id[a]] = map[a]; 730 772 } 731 773 return *this; 732 }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 before749 /// calling \ref run(), the lower bounds will be set to zero750 /// 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 to754 /// \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 the760 /// \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 774 } 770 775 … … 780 785 /// 781 786 /// \return <tt>(*this)</tt> 782 template<typename COST> 783 NetworkSimplex& costMap(const COST& map) { 784 delete _pcost; 785 _pcost = new CostArcMap(_graph); 787 template<typename CostMap> 788 NetworkSimplex& costMap(const CostMap& map) { 786 789 for (ArcIt a(_graph); a != INVALID; ++a) { 787 (*_pcost)[a] = map[a];790 _cost[_arc_id[a]] = map[a]; 788 791 } 789 792 return *this; … … 798 801 /// 799 802 /// \param map A node map storing the supply values. 800 /// Its \c Value type must be convertible to the \c Flowtype803 /// Its \c Value type must be convertible to the \c Value type 801 804 /// of the algorithm. 802 805 /// 803 806 /// \return <tt>(*this)</tt> 804 template<typename SUP> 805 NetworkSimplex& supplyMap(const SUP& map) { 806 delete _psupply; 807 _pstsup = false; 808 _psupply = new FlowNodeMap(_graph); 807 template<typename SupplyMap> 808 NetworkSimplex& supplyMap(const SupplyMap& map) { 809 809 for (NodeIt n(_graph); n != INVALID; ++n) { 810 (*_psupply)[n] = map[n];810 _supply[_node_id[n]] = map[n]; 811 811 } 812 812 return *this; … … 821 821 /// (It makes sense only if non-zero lower bounds are given.) 822 822 /// 823 /// Using this function has the same effect as using \ref supplyMap() 824 /// with such a map in which \c k is assigned to \c s, \c -k is 825 /// assigned to \c t and all other nodes have zero supply value. 826 /// 823 827 /// \param s The source node. 824 828 /// \param t The target node. … … 827 831 /// 828 832 /// \return <tt>(*this)</tt> 829 NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) { 830 delete _psupply; 831 _psupply = NULL; 832 _pstsup = true; 833 _psource = s; 834 _ptarget = t; 835 _pstflow = k; 833 NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) { 834 for (int i = 0; i != _node_num; ++i) { 835 _supply[i] = 0; 836 } 837 _supply[_node_id[s]] = k; 838 _supply[_node_id[t]] = -k; 836 839 return *this; 837 840 } 838 841 839 /// \brief Set the problem type.840 /// 841 /// This function sets the problem type for the algorithm.842 /// If it is not used before calling \ref run(), the \ref GEQ problem842 /// \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 843 846 /// type will be used. 844 847 /// 845 /// For more information see \ref ProblemType.848 /// For more information see \ref SupplyType. 846 849 /// 847 850 /// \return <tt>(*this)</tt> 848 NetworkSimplex& problemType(ProblemType problem_type) {849 _ ptype = problem_type;851 NetworkSimplex& supplyType(SupplyType supply_type) { 852 _stype = supply_type; 850 853 return *this; 851 854 } 852 855 853 /// \brief Set the flow map.854 ///855 /// This function sets the flow map.856 /// If it is not used before calling \ref run(), an instance will857 /// be allocated automatically. The destructor deallocates this858 /// 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 storing873 /// the dual solution.874 /// If it is not used before calling \ref run(), an instance will875 /// be allocated automatically. The destructor deallocates this876 /// 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 856 /// @} 889 857 … … 897 865 /// This function runs the algorithm. 898 866 /// The paramters can be specified using functions \ref lowerMap(), 899 /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(), 900 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 901 /// \ref problemType(), \ref flowMap() and \ref potentialMap(). 867 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 868 /// \ref supplyType(). 902 869 /// For example, 903 870 /// \code 904 871 /// NetworkSimplex<ListDigraph> ns(graph); 905 /// ns. boundMaps(lower,upper).costMap(cost)872 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 906 873 /// .supplyMap(sup).run(); 907 874 /// \endcode … … 911 878 /// \ref reset() is called, thus only the modified parameters 912 879 /// have to be set again. See \ref reset() for examples. 880 /// However the underlying digraph must not be modified after this 881 /// class have been constructed, since it copies and extends the graph. 913 882 /// 914 883 /// \param pivot_rule The pivot rule that will be used during the 915 884 /// algorithm. For more information see \ref PivotRule. 916 885 /// 917 /// \return \c true if a feasible flow can be found. 918 bool run(PivotRule pivot_rule = BLOCK_SEARCH) { 919 return init() && start(pivot_rule); 886 /// \return \c INFEASIBLE if no feasible flow exists, 887 /// \n \c OPTIMAL if the problem has optimal solution 888 /// (i.e. it is feasible and bounded), and the algorithm has found 889 /// optimal flow and node potentials (primal and dual solutions), 890 /// \n \c UNBOUNDED if the objective function of the problem is 891 /// unbounded, i.e. there is a directed cycle having negative total 892 /// cost and infinite upper bound. 893 /// 894 /// \see ProblemType, PivotRule 895 ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) { 896 if (!init()) return INFEASIBLE; 897 return start(pivot_rule); 920 898 } 921 899 … … 924 902 /// This function resets all the paramaters that have been given 925 903 /// before using functions \ref lowerMap(), \ref upperMap(), 926 /// \ref capacityMap(), \ref boundMaps(), \ref costMap(), 927 /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 928 /// \ref flowMap() and \ref potentialMap(). 904 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(). 929 905 /// 930 906 /// It is useful for multiple run() calls. If this function is not 931 907 /// used, all the parameters given before are kept for the next 932 908 /// \ref run() call. 909 /// However the underlying digraph must not be modified after this 910 /// class have been constructed, since it copies and extends the graph. 933 911 /// 934 912 /// For example, … … 937 915 /// 938 916 /// // First run 939 /// ns.lowerMap(lower). capacityMap(cap).costMap(cost)917 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 940 918 /// .supplyMap(sup).run(); 941 919 /// … … 948 926 /// // (the lower bounds will be set to zero on all arcs) 949 927 /// ns.reset(); 950 /// ns. capacityMap(cap).costMap(cost)928 /// ns.upperMap(capacity).costMap(cost) 951 929 /// .supplyMap(sup).run(); 952 930 /// \endcode … … 954 932 /// \return <tt>(*this)</tt> 955 933 NetworkSimplex& reset() { 956 delete _plower; 957 delete _pupper; 958 delete _pcost; 959 delete _psupply; 960 _plower = NULL; 961 _pupper = NULL; 962 _pcost = NULL; 963 _psupply = NULL; 964 _pstsup = false; 965 _ptype = GEQ; 966 if (_local_flow) delete _flow_map; 967 if (_local_potential) delete _potential_map; 968 _flow_map = NULL; 969 _potential_map = NULL; 970 _local_flow = false; 971 _local_potential = false; 972 934 for (int i = 0; i != _node_num; ++i) { 935 _supply[i] = 0; 936 } 937 for (int i = 0; i != _arc_num; ++i) { 938 _lower[i] = 0; 939 _upper[i] = INF; 940 _cost[i] = 1; 941 } 942 _have_lower = false; 943 _stype = GEQ; 973 944 return *this; 974 945 } … … 986 957 /// 987 958 /// This function returns the total cost of the found flow. 988 /// The complexity of the functionis O(e).959 /// Its complexity is O(e). 989 960 /// 990 961 /// \note The return type of the function can be specified as a … … 998 969 /// 999 970 /// \pre \ref run() must be called before using this function. 1000 template <typename Num> 1001 Num totalCost() const { 1002 Num c = 0; 1003 if (_pcost) { 1004 for (ArcIt e(_graph); e != INVALID; ++e) 1005 c += (*_flow_map)[e] * (*_pcost)[e]; 1006 } else { 1007 for (ArcIt e(_graph); e != INVALID; ++e) 1008 c += (*_flow_map)[e]; 971 template <typename Number> 972 Number totalCost() const { 973 Number c = 0; 974 for (ArcIt a(_graph); a != INVALID; ++a) { 975 int i = _arc_id[a]; 976 c += Number(_flow[i]) * Number(_cost[i]); 1009 977 } 1010 978 return c; … … 1022 990 /// 1023 991 /// \pre \ref run() must be called before using this function. 1024 Flow flow(const Arc& a) const { 1025 return (*_flow_map)[a]; 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. 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. 1032 1001 /// 1033 1002 /// \pre \ref run() must be called before using this function. 1034 const FlowMap& flowMap() const { 1035 return *_flow_map; 1003 template <typename FlowMap> 1004 void flowMap(FlowMap &map) const { 1005 for (ArcIt a(_graph); a != INVALID; ++a) { 1006 map.set(a, _flow[_arc_id[a]]); 1007 } 1036 1008 } 1037 1009 … … 1043 1015 /// \pre \ref run() must be called before using this function. 1044 1016 Cost potential(const Node& n) const { 1045 return (*_potential_map)[n];1046 } 1047 1048 /// \brief Return a const reference to the potential map1049 /// (the dual solution).1050 /// 1051 /// This function returns a const reference to a node map storing1052 /// the found potentials, which form the dual solution ofthe1053 /// \ ref min_cost_flow "minimum cost flow" problem.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. 1054 1026 /// 1055 1027 /// \pre \ref run() must be called before using this function. 1056 const PotentialMap& potentialMap() const { 1057 return *_potential_map; 1028 template <typename PotentialMap> 1029 void potentialMap(PotentialMap &map) const { 1030 for (NodeIt n(_graph); n != INVALID; ++n) { 1031 map.set(n, _pi[_node_id[n]]); 1032 } 1058 1033 } 1059 1034 … … 1064 1039 // Initialize internal data structures 1065 1040 bool init() { 1066 // Initialize result maps1067 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 vectors1077 _node_num = countNodes(_graph);1078 _arc_num = countArcs(_graph);1079 int all_node_num = _node_num + 1;1080 int all_arc_num = _arc_num + _node_num;1081 1041 if (_node_num == 0) return false; 1082 1042 1083 _arc_ref.resize(_arc_num); 1084 _source.resize(all_arc_num); 1085 _target.resize(all_arc_num); 1086 1087 _cap.resize(all_arc_num); 1088 _cost.resize(all_arc_num); 1089 _supply.resize(all_node_num); 1090 _flow.resize(all_arc_num); 1091 _pi.resize(all_node_num); 1092 1093 _parent.resize(all_node_num); 1094 _pred.resize(all_node_num); 1095 _forward.resize(all_node_num); 1096 _thread.resize(all_node_num); 1097 _rev_thread.resize(all_node_num); 1098 _succ_num.resize(all_node_num); 1099 _last_succ.resize(all_node_num); 1100 _state.resize(all_arc_num); 1101 1102 // Initialize node related data 1103 bool valid_supply = true; 1104 Flow sum_supply = 0; 1105 if (!_pstsup && !_psupply) { 1106 _pstsup = true; 1107 _psource = _ptarget = NodeIt(_graph); 1108 _pstflow = 0; 1109 } 1110 if (_psupply) { 1111 int i = 0; 1112 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 1113 _node_id[n] = i; 1114 _supply[i] = (*_psupply)[n]; 1115 sum_supply += _supply[i]; 1116 } 1117 valid_supply = (_ptype == GEQ && sum_supply <= 0) || 1118 (_ptype == LEQ && sum_supply >= 0); 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) { 1053 for (int i = 0; i != _arc_num; ++i) { 1054 Value c = _lower[i]; 1055 if (c >= 0) { 1056 _cap[i] = _upper[i] < INF ? _upper[i] - c : INF; 1057 } else { 1058 _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF; 1059 } 1060 _supply[_source[i]] -= c; 1061 _supply[_target[i]] += c; 1062 } 1119 1063 } 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(); 1064 for (int i = 0; i != _arc_num; ++i) { 1065 _cap[i] = _upper[i]; 1066 } 1067 } 1135 1068 1136 1069 // Initialize artifical cost 1137 Cost art_cost;1070 Cost ART_COST; 1138 1071 if (std::numeric_limits<Cost>::is_exact) { 1139 art_cost= std::numeric_limits<Cost>::max() / 4 + 1;1072 ART_COST = std::numeric_limits<Cost>::max() / 4 + 1; 1140 1073 } else { 1141 art_cost= std::numeric_limits<Cost>::min();1074 ART_COST = std::numeric_limits<Cost>::min(); 1142 1075 for (int i = 0; i != _arc_num; ++i) { 1143 if (_cost[i] > art_cost) art_cost = _cost[i]; 1144 } 1145 art_cost = (art_cost + 1) * _node_num; 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(); 1169 } else { 1170 Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap> 1171 circ(_graph, *_plower, inf_arc_map, *csup); 1172 circ_result = circ.run(); 1173 } 1174 } else { 1175 if (_pupper) { 1176 Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap> 1177 circ(_graph, zero_arc_map, *_pupper, *csup); 1178 circ_result = circ.run(); 1179 } else { 1180 Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap> 1181 circ(_graph, zero_arc_map, inf_arc_map, *csup); 1182 circ_result = circ.run(); 1183 } 1184 } 1185 } else { 1186 // LEQ problem type 1187 typedef ReverseDigraph<const GR> RevGraph; 1188 typedef NegMap<FlowNodeMap> NegNodeMap; 1189 RevGraph rgraph(_graph); 1190 NegNodeMap neg_csup(*csup); 1191 if (_plower) { 1192 if (_pupper) { 1193 Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap> 1194 circ(rgraph, *_plower, *_pupper, neg_csup); 1195 circ_result = circ.run(); 1196 } else { 1197 Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap> 1198 circ(rgraph, *_plower, inf_arc_map, neg_csup); 1199 circ_result = circ.run(); 1200 } 1201 } else { 1202 if (_pupper) { 1203 Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap> 1204 circ(rgraph, zero_arc_map, *_pupper, neg_csup); 1205 circ_result = circ.run(); 1206 } else { 1207 Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap> 1208 circ(rgraph, zero_arc_map, inf_arc_map, neg_csup); 1209 circ_result = circ.run(); 1210 } 1211 } 1212 } 1213 if (local_csup) delete csup; 1214 if (!circ_result) return false; 1215 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 1216 1087 // Set data for the artificial root node 1217 1088 _root = _node_num; … … 1220 1091 _thread[_root] = 0; 1221 1092 _rev_thread[0] = _root; 1222 _succ_num[_root] = all_node_num;1093 _succ_num[_root] = _node_num + 1; 1223 1094 _last_succ[_root] = _root - 1; 1224 _supply[_root] = -sum_supply; 1225 if (sum_supply < 0) { 1226 _pi[_root] = -art_cost; 1227 } else { 1228 _pi[_root] = art_cost; 1229 } 1230 1231 // Store the arcs in a mixed order 1232 int k = std::max(int(std::sqrt(double(_arc_num))), 10); 1233 int i = 0; 1234 for (ArcIt e(_graph); e != INVALID; ++e) { 1235 _arc_ref[i] = e; 1236 if ((i += k) >= _arc_num) i = (i % k) + 1; 1237 } 1238 1239 // Initialize arc maps 1240 if (_pupper && _pcost) { 1241 for (int i = 0; i != _arc_num; ++i) { 1242 Arc e = _arc_ref[i]; 1243 _source[i] = _node_id[_graph.source(e)]; 1244 _target[i] = _node_id[_graph.target(e)]; 1245 _cap[i] = (*_pupper)[e]; 1246 _cost[i] = (*_pcost)[e]; 1247 _flow[i] = 0; 1248 _state[i] = STATE_LOWER; 1249 } 1250 } else { 1251 for (int i = 0; i != _arc_num; ++i) { 1252 Arc e = _arc_ref[i]; 1253 _source[i] = _node_id[_graph.source(e)]; 1254 _target[i] = _node_id[_graph.target(e)]; 1255 _flow[i] = 0; 1256 _state[i] = STATE_LOWER; 1257 } 1258 if (_pupper) { 1259 for (int i = 0; i != _arc_num; ++i) 1260 _cap[i] = (*_pupper)[_arc_ref[i]]; 1261 } else { 1262 for (int i = 0; i != _arc_num; ++i) 1263 _cap[i] = inf_cap; 1264 } 1265 if (_pcost) { 1266 for (int i = 0; i != _arc_num; ++i) 1267 _cost[i] = (*_pcost)[_arc_ref[i]]; 1268 } else { 1269 for (int i = 0; i != _arc_num; ++i) 1270 _cost[i] = 1; 1271 } 1272 } 1273 1274 // Remove non-zero lower bounds 1275 if (_plower) { 1276 for (int i = 0; i != _arc_num; ++i) { 1277 Flow c = (*_plower)[_arc_ref[i]]; 1278 if (c != 0) { 1279 _cap[i] -= c; 1280 _supply[_source[i]] -= c; 1281 _supply[_target[i]] += c; 1282 } 1283 } 1284 } 1095 _supply[_root] = -_sum_supply; 1096 _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST; 1285 1097 1286 1098 // Add artificial arcs and initialize the spanning tree data structure 1287 1099 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1100 _parent[u] = _root; 1101 _pred[u] = e; 1288 1102 _thread[u] = u + 1; 1289 1103 _rev_thread[u + 1] = u; 1290 1104 _succ_num[u] = 1; 1291 1105 _last_succ[u] = u; 1292 _parent[u] = _root; 1293 _pred[u] = e; 1294 _cost[e] = art_cost; 1295 _cap[e] = inf_cap; 1106 _cost[e] = ART_COST; 1107 _cap[e] = INF; 1296 1108 _state[e] = STATE_TREE; 1297 if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {1109 if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) { 1298 1110 _flow[e] = _supply[u]; 1299 1111 _forward[u] = true; 1300 _pi[u] = - art_cost+ _pi[_root];1112 _pi[u] = -ART_COST + _pi[_root]; 1301 1113 } else { 1302 1114 _flow[e] = -_supply[u]; 1303 1115 _forward[u] = false; 1304 _pi[u] = art_cost+ _pi[_root];1116 _pi[u] = ART_COST + _pi[_root]; 1305 1117 } 1306 1118 } … … 1337 1149 delta = _cap[in_arc]; 1338 1150 int result = 0; 1339 Flowd;1151 Value d; 1340 1152 int e; 1341 1153 … … 1343 1155 for (int u = first; u != join; u = _parent[u]) { 1344 1156 e = _pred[u]; 1345 d = _forward[u] ? _flow[e] : _cap[e] - _flow[e]; 1157 d = _forward[u] ? 1158 _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]); 1346 1159 if (d < delta) { 1347 1160 delta = d; … … 1353 1166 for (int u = second; u != join; u = _parent[u]) { 1354 1167 e = _pred[u]; 1355 d = _forward[u] ? _cap[e] - _flow[e] : _flow[e]; 1168 d = _forward[u] ? 1169 (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e]; 1356 1170 if (d <= delta) { 1357 1171 delta = d; … … 1375 1189 // Augment along the cycle 1376 1190 if (delta > 0) { 1377 Flowval = _state[in_arc] * delta;1191 Value val = _state[in_arc] * delta; 1378 1192 _flow[in_arc] += val; 1379 1193 for (int u = _source[in_arc]; u != join; u = _parent[u]) { … … 1527 1341 1528 1342 // Execute the algorithm 1529 boolstart(PivotRule pivot_rule) {1343 ProblemType start(PivotRule pivot_rule) { 1530 1344 // Select the pivot rule implementation 1531 1345 switch (pivot_rule) { … … 1541 1355 return start<AlteringListPivotRule>(); 1542 1356 } 1543 return false;1357 return INFEASIBLE; // avoid warning 1544 1358 } 1545 1359 1546 1360 template <typename PivotRuleImpl> 1547 boolstart() {1361 ProblemType start() { 1548 1362 PivotRuleImpl pivot(*this); 1549 1363 … … 1552 1366 findJoinNode(); 1553 1367 bool change = findLeavingArc(); 1368 if (delta >= INF) return UNBOUNDED; 1554 1369 changeFlow(change); 1555 1370 if (change) { … … 1558 1373 } 1559 1374 } 1560 1561 // Copy flow values to _flow_map 1562 if (_plower) { 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) { 1563 1395 for (int i = 0; i != _arc_num; ++i) { 1564 Arc e = _arc_ref[i]; 1565 _flow_map->set(e, (*_plower)[e] + _flow[i]); 1566 } 1567 } else { 1568 for (int i = 0; i != _arc_num; ++i) { 1569 _flow_map->set(_arc_ref[i], _flow[i]); 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; 1396 Value c = _lower[i]; 1397 if (c != 0) { 1398 _flow[i] += c; 1399 _supply[_source[i]] += c; 1400 _supply[_target[i]] -= c; 1401 } 1402 } 1403 } 1404 1405 return OPTIMAL; 1578 1406 } 1579 1407 -
lemon/preflow.h
r658 r688 47 47 48 48 /// \brief The type of the flow values. 49 typedef typename CapacityMap::Value Flow;49 typedef typename CapacityMap::Value Value; 50 50 51 51 /// \brief The type of the map that stores the flow values. … … 53 53 /// The type of the map that stores the flow values. 54 54 /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept. 55 typedef typename Digraph::template ArcMap< Flow> FlowMap;55 typedef typename Digraph::template ArcMap<Value> FlowMap; 56 56 57 57 /// \brief Instantiates a FlowMap. … … 85 85 /// 86 86 /// The tolerance used by the algorithm to handle inexact computation. 87 typedef lemon::Tolerance< Flow> Tolerance;87 typedef lemon::Tolerance<Value> Tolerance; 88 88 89 89 }; … … 126 126 typedef typename Traits::CapacityMap CapacityMap; 127 127 ///The type of the flow values. 128 typedef typename Traits:: Flow Flow;128 typedef typename Traits::Value Value; 129 129 130 130 ///The type of the flow map. … … 152 152 bool _local_level; 153 153 154 typedef typename Digraph::template NodeMap< Flow> ExcessMap;154 typedef typename Digraph::template NodeMap<Value> ExcessMap; 155 155 ExcessMap* _excess; 156 156 … … 471 471 472 472 for (NodeIt n(_graph); n != INVALID; ++n) { 473 Flowexcess = 0;473 Value 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 Flowrem = (*_capacity)[e] - (*_flow)[e];522 Value 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 Flowrem = (*_flow)[e];534 Value 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 Flowexcess = (*_excess)[n];567 Value excess = (*_excess)[n]; 568 568 int new_level = _level->maxLevel(); 569 569 570 570 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 571 Flowrem = (*_capacity)[e] - (*_flow)[e];571 Value rem = (*_capacity)[e] - (*_flow)[e]; 572 572 if (!_tolerance.positive(rem)) continue; 573 573 Node v = _graph.target(e); … … 592 592 593 593 for (InArcIt e(_graph, n); e != INVALID; ++e) { 594 Flowrem = (*_flow)[e];594 Value 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 Flowexcess = (*_excess)[n];640 Value excess = (*_excess)[n]; 641 641 int new_level = _level->maxLevel(); 642 642 643 643 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 644 Flowrem = (*_capacity)[e] - (*_flow)[e];644 Value rem = (*_capacity)[e] - (*_flow)[e]; 645 645 if (!_tolerance.positive(rem)) continue; 646 646 Node v = _graph.target(e); … … 665 665 666 666 for (InArcIt e(_graph, n); e != INVALID; ++e) { 667 Flowrem = (*_flow)[e];667 Value 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 Flowexcess = (*_excess)[n];781 Value excess = (*_excess)[n]; 782 782 int level = _level->highestActiveLevel(); 783 783 int new_level = _level->maxLevel(); 784 784 785 785 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 786 Flowrem = (*_capacity)[e] - (*_flow)[e];786 Value rem = (*_capacity)[e] - (*_flow)[e]; 787 787 if (!_tolerance.positive(rem)) continue; 788 788 Node v = _graph.target(e); … … 807 807 808 808 for (InArcIt e(_graph, n); e != INVALID; ++e) { 809 Flowrem = (*_flow)[e];809 Value 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 FlowflowValue() const {900 Value flowValue() const { 901 901 return (*_excess)[_target]; 902 902 } 903 903 904 /// \brief Returns the flow on the given arc.905 /// 906 /// Returns the flow on the given arc. This method can904 /// \brief Returns the flow value on the given arc. 905 /// 906 /// Returns the flow value on the given arc. This method can 907 907 /// be called after the second phase of the algorithm. 908 908 /// 909 909 /// \pre Either \ref run() or \ref init() must be called before 910 910 /// using this function. 911 Flowflow(const Arc& arc) const {911 Value flow(const Arc& arc) const { 912 912 return (*_flow)[arc]; 913 913 } -
test/lp_test.cc
r674 r678 22 22 #include <lemon/tolerance.h> 23 23 24 #ifdef HAVE_CONFIG_H25 24 #include <lemon/config.h> 26 #endif27 25 28 26 #ifdef LEMON_HAVE_GLPK -
test/min_cost_flow_test.cc
r662 r689 19 19 #include <iostream> 20 20 #include <fstream> 21 #include <limits> 21 22 22 23 #include <lemon/list_graph.h> … … 34 35 char test_lgf[] = 35 36 "@nodes\n" 36 "label sup1 sup2 sup3 sup4 sup5 \n"37 " 1 20 27 0 20 30\n"38 " 2 -4 0 0 -8 -3\n"39 " 3 0 0 0 0 0 \n"40 " 4 0 0 0 0 0 \n"41 " 5 9 0 0 6 11\n"42 " 6 -6 0 0 -5 -6\n"43 " 7 0 0 0 0 0 \n"44 " 8 0 0 0 0 3\n"45 " 9 3 0 0 0 0 \n"46 " 10 -2 0 0 -7 -2\n"47 " 11 0 0 0 -10 0\n"48 " 12 -20 -27 0 -30 - 20\n"49 "\n" 37 "label sup1 sup2 sup3 sup4 sup5 sup6\n" 38 " 1 20 27 0 30 20 30\n" 39 " 2 -4 0 0 0 -8 -3\n" 40 " 3 0 0 0 0 0 0\n" 41 " 4 0 0 0 0 0 0\n" 42 " 5 9 0 0 0 6 11\n" 43 " 6 -6 0 0 0 -5 -6\n" 44 " 7 0 0 0 0 0 0\n" 45 " 8 0 0 0 0 0 3\n" 46 " 9 3 0 0 0 0 0\n" 47 " 10 -2 0 0 0 -7 -2\n" 48 " 11 0 0 0 0 -10 0\n" 49 " 12 -20 -27 0 -30 -30 -20\n" 50 "\n" 50 51 "@arcs\n" 51 " cost cap low1 low2 \n"52 " 1 2 70 11 0 8 \n"53 " 1 3 150 3 0 1 \n"54 " 1 4 80 15 0 2 \n"55 " 2 8 80 12 0 0 \n"56 " 3 5 140 5 0 3 \n"57 " 4 6 60 10 0 1 \n"58 " 4 7 80 2 0 0 \n"59 " 4 8 110 3 0 0 \n"60 " 5 7 60 14 0 0 \n"61 " 5 11 120 12 0 0 \n"62 " 6 3 0 3 0 0 \n"63 " 6 9 140 4 0 0 \n"64 " 6 10 90 8 0 0 \n"65 " 7 1 30 5 0 0 \n"66 " 8 12 60 16 0 4 \n"67 " 9 12 50 6 0 0 \n"68 "10 12 70 13 0 5 \n"69 "10 2 100 7 0 0 \n"70 "10 7 60 10 0 0 \n"71 "11 10 20 14 0 6 \n"72 "12 11 30 10 0 0 \n"52 " cost cap low1 low2 low3\n" 53 " 1 2 70 11 0 8 8\n" 54 " 1 3 150 3 0 1 0\n" 55 " 1 4 80 15 0 2 2\n" 56 " 2 8 80 12 0 0 0\n" 57 " 3 5 140 5 0 3 1\n" 58 " 4 6 60 10 0 1 0\n" 59 " 4 7 80 2 0 0 0\n" 60 " 4 8 110 3 0 0 0\n" 61 " 5 7 60 14 0 0 0\n" 62 " 5 11 120 12 0 0 0\n" 63 " 6 3 0 3 0 0 0\n" 64 " 6 9 140 4 0 0 0\n" 65 " 6 10 90 8 0 0 0\n" 66 " 7 1 30 5 0 0 -5\n" 67 " 8 12 60 16 0 4 3\n" 68 " 9 12 50 6 0 0 0\n" 69 "10 12 70 13 0 5 2\n" 70 "10 2 100 7 0 0 0\n" 71 "10 7 60 10 0 0 -3\n" 72 "11 10 20 14 0 6 -20\n" 73 "12 11 30 10 0 0 -10\n" 73 74 "\n" 74 75 "@attributes\n" … … 77 78 78 79 79 enum ProblemType {80 enum SupplyType { 80 81 EQ, 81 82 GEQ, … … 84 85 85 86 // Check the interface of an MCF algorithm 86 template <typename GR, typename Flow, typename Cost>87 template <typename GR, typename Value, typename Cost> 87 88 class McfClassConcept 88 89 { … … 95 96 96 97 MCF mcf(g); 98 const MCF& const_mcf = mcf; 97 99 98 100 b = mcf.reset() 99 101 .lowerMap(lower) 100 102 .upperMap(upper) 101 .capacityMap(upper)102 .boundMaps(lower, upper)103 103 .costMap(cost) 104 104 .supplyMap(sup) 105 105 .stSupply(n, n, k) 106 .flowMap(flow)107 .potentialMap(pot)108 106 .run(); 109 110 const MCF& const_mcf = mcf; 111 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>(); 107 108 c = const_mcf.totalCost(); 109 x = const_mcf.template totalCost<double>(); 117 110 v = const_mcf.flow(a); 118 v = const_mcf.potential(n); 119 120 ignore_unused_variable_warning(fm); 121 ignore_unused_variable_warning(pm); 122 ignore_unused_variable_warning(x); 111 c = const_mcf.potential(n); 112 const_mcf.flowMap(fm); 113 const_mcf.potentialMap(pm); 123 114 } 124 115 125 116 typedef typename GR::Node Node; 126 117 typedef typename GR::Arc Arc; 127 typedef concepts::ReadMap<Node, Flow> NM;128 typedef concepts::ReadMap<Arc, Flow> FAM;118 typedef concepts::ReadMap<Node, Value> NM; 119 typedef concepts::ReadMap<Arc, Value> VAM; 129 120 typedef concepts::ReadMap<Arc, Cost> CAM; 121 typedef concepts::WriteMap<Arc, Value> FlowMap; 122 typedef concepts::WriteMap<Node, Cost> PotMap; 130 123 131 124 const GR &g; 132 const FAM &lower;133 const FAM &upper;125 const VAM &lower; 126 const VAM &upper; 134 127 const CAM &cost; 135 128 const NM ⊃ 136 129 const Node &n; 137 130 const Arc &a; 138 const Flow &k; 139 Flow v; 131 const Value &k; 132 FlowMap fm; 133 PotMap pm; 140 134 bool b; 141 142 typename MCF:: FlowMap &flow;143 typename MCF:: PotentialMap &pot;135 double x; 136 typename MCF::Value v; 137 typename MCF::Cost c; 144 138 }; 145 139 … … 152 146 bool checkFlow( const GR& gr, const LM& lower, const UM& upper, 153 147 const SM& supply, const FM& flow, 154 ProblemType type = EQ )148 SupplyType type = EQ ) 155 149 { 156 150 TEMPLATE_DIGRAPH_TYPEDEFS(GR); … … 209 203 template < typename MCF, typename GR, 210 204 typename LM, typename UM, 211 typename CM, typename SM > 212 void checkMcf( const MCF& mcf, bool mcf_result, 205 typename CM, typename SM, 206 typename PT > 207 void checkMcf( const MCF& mcf, PT mcf_result, 213 208 const GR& gr, const LM& lower, const UM& upper, 214 209 const CM& cost, const SM& supply, 215 bool result, typename CM::Value total,210 PT result, bool optimal, typename CM::Value total, 216 211 const std::string &test_id = "", 217 ProblemType type = EQ )212 SupplyType type = EQ ) 218 213 { 219 214 check(mcf_result == result, "Wrong result " + test_id); 220 if (result) { 221 check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type), 215 if (optimal) { 216 typename GR::template ArcMap<typename SM::Value> flow(gr); 217 typename GR::template NodeMap<typename CM::Value> pi(gr); 218 mcf.flowMap(flow); 219 mcf.potentialMap(pi); 220 check(checkFlow(gr, lower, upper, supply, flow, type), 222 221 "The flow is not feasible " + test_id); 223 222 check(mcf.totalCost() == total, "The flow is not optimal " + test_id); 224 check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(), 225 mcf.potentialMap()), 223 check(checkPotential(gr, lower, upper, cost, supply, flow, pi), 226 224 "Wrong potentials " + test_id); 227 225 } … … 232 230 // Check the interfaces 233 231 { 234 typedef int Flow;235 typedef int Cost;236 232 typedef concepts::Digraph GR; 237 checkConcept< McfClassConcept<GR, Flow, Cost>, 238 NetworkSimplex<GR, Flow, Cost> >(); 233 checkConcept< McfClassConcept<GR, int, int>, 234 NetworkSimplex<GR> >(); 235 checkConcept< McfClassConcept<GR, double, double>, 236 NetworkSimplex<GR, double> >(); 237 checkConcept< McfClassConcept<GR, int, double>, 238 NetworkSimplex<GR, int, double> >(); 239 239 } 240 240 … … 245 245 // Read the test digraph 246 246 Digraph gr; 247 Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);248 Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr) ;247 Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr); 248 Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr); 249 249 ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max()); 250 250 Node v, w; … … 256 256 .arcMap("low1", l1) 257 257 .arcMap("low2", l2) 258 .arcMap("low3", l3) 258 259 .nodeMap("sup1", s1) 259 260 .nodeMap("sup2", s2) … … 261 262 .nodeMap("sup4", s4) 262 263 .nodeMap("sup5", s5) 264 .nodeMap("sup6", s6) 263 265 .node("source", v) 264 266 .node("target", w) 265 267 .run(); 268 269 // Build a test digraph for testing negative costs 270 Digraph ngr; 271 Node n1 = ngr.addNode(); 272 Node n2 = ngr.addNode(); 273 Node n3 = ngr.addNode(); 274 Node n4 = ngr.addNode(); 275 Node n5 = ngr.addNode(); 276 Node n6 = ngr.addNode(); 277 Node n7 = ngr.addNode(); 278 279 Arc a1 = ngr.addArc(n1, n2); 280 Arc a2 = ngr.addArc(n1, n3); 281 Arc a3 = ngr.addArc(n2, n4); 282 Arc a4 = ngr.addArc(n3, n4); 283 Arc a5 = ngr.addArc(n3, n2); 284 Arc a6 = ngr.addArc(n5, n3); 285 Arc a7 = ngr.addArc(n5, n6); 286 Arc a8 = ngr.addArc(n6, n7); 287 Arc a9 = ngr.addArc(n7, n5); 288 289 Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0); 290 ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000); 291 Digraph::NodeMap<int> ns(ngr, 0); 292 293 nl2[a7] = 1000; 294 nl2[a8] = -1000; 295 296 ns[n1] = 100; 297 ns[n4] = -100; 298 299 nc[a1] = 100; 300 nc[a2] = 30; 301 nc[a3] = 20; 302 nc[a4] = 80; 303 nc[a5] = 50; 304 nc[a6] = 10; 305 nc[a7] = 80; 306 nc[a8] = 30; 307 nc[a9] = -120; 266 308 267 309 // A. Test NetworkSimplex with the default pivot rule … … 272 314 mcf.upperMap(u).costMap(c); 273 315 checkMcf(mcf, mcf.supplyMap(s1).run(), 274 gr, l1, u, c, s1, true,5240, "#A1");316 gr, l1, u, c, s1, mcf.OPTIMAL, true, 5240, "#A1"); 275 317 checkMcf(mcf, mcf.stSupply(v, w, 27).run(), 276 gr, l1, u, c, s2, true,7620, "#A2");318 gr, l1, u, c, s2, mcf.OPTIMAL, true, 7620, "#A2"); 277 319 mcf.lowerMap(l2); 278 320 checkMcf(mcf, mcf.supplyMap(s1).run(), 279 gr, l2, u, c, s1, true,5970, "#A3");321 gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#A3"); 280 322 checkMcf(mcf, mcf.stSupply(v, w, 27).run(), 281 gr, l2, u, c, s2, true,8010, "#A4");323 gr, l2, u, c, s2, mcf.OPTIMAL, true, 8010, "#A4"); 282 324 mcf.reset(); 283 325 checkMcf(mcf, mcf.supplyMap(s1).run(), 284 gr, l1, cu, cc, s1, true,74, "#A5");326 gr, l1, cu, cc, s1, mcf.OPTIMAL, true, 74, "#A5"); 285 327 checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(), 286 gr, l2, cu, cc, s2, true,94, "#A6");328 gr, l2, cu, cc, s2, mcf.OPTIMAL, true, 94, "#A6"); 287 329 mcf.reset(); 288 330 checkMcf(mcf, mcf.run(), 289 gr, l1, cu, cc, s3, true, 0, "#A7"); 290 checkMcf(mcf, mcf.boundMaps(l2, u).run(), 291 gr, l2, u, cc, s3, false, 0, "#A8"); 331 gr, l1, cu, cc, s3, mcf.OPTIMAL, true, 0, "#A7"); 332 checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(), 333 gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8"); 334 mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4); 335 checkMcf(mcf, mcf.run(), 336 gr, l3, u, c, s4, mcf.OPTIMAL, true, 6360, "#A9"); 292 337 293 338 // Check the GEQ form 294 mcf.reset().upperMap(u).costMap(c).supplyMap(s 4);295 checkMcf(mcf, mcf.run(), 296 gr, l1, u, c, s 4, true, 3530, "#A9", GEQ);297 mcf. problemType(mcf.GEQ);339 mcf.reset().upperMap(u).costMap(c).supplyMap(s5); 340 checkMcf(mcf, mcf.run(), 341 gr, l1, u, c, s5, mcf.OPTIMAL, true, 3530, "#A10", GEQ); 342 mcf.supplyType(mcf.GEQ); 298 343 checkMcf(mcf, mcf.lowerMap(l2).run(), 299 gr, l2, u, c, s 4, true, 4540, "#A10", GEQ);300 mcf. problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);301 checkMcf(mcf, mcf.run(), 302 gr, l2, u, c, s 5, false, 0, "#A11", GEQ);344 gr, l2, u, c, s5, mcf.OPTIMAL, true, 4540, "#A11", GEQ); 345 mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6); 346 checkMcf(mcf, mcf.run(), 347 gr, l2, u, c, s6, mcf.INFEASIBLE, false, 0, "#A12", GEQ); 303 348 304 349 // Check the LEQ form 305 mcf.reset(). problemType(mcf.LEQ);306 mcf.upperMap(u).costMap(c).supplyMap(s 5);307 checkMcf(mcf, mcf.run(), 308 gr, l1, u, c, s 5, true, 5080, "#A12", LEQ);350 mcf.reset().supplyType(mcf.LEQ); 351 mcf.upperMap(u).costMap(c).supplyMap(s6); 352 checkMcf(mcf, mcf.run(), 353 gr, l1, u, c, s6, mcf.OPTIMAL, true, 5080, "#A13", LEQ); 309 354 checkMcf(mcf, mcf.lowerMap(l2).run(), 310 gr, l2, u, c, s5, true, 5930, "#A13", LEQ); 311 mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4); 312 checkMcf(mcf, mcf.run(), 313 gr, l2, u, c, s4, false, 0, "#A14", LEQ); 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"); 314 370 } 315 371 … … 317 373 { 318 374 NetworkSimplex<Digraph> mcf(gr); 319 mcf.supplyMap(s1).costMap(c). capacityMap(u).lowerMap(l2);375 mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2); 320 376 321 377 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE), 322 gr, l2, u, c, s1, true,5970, "#B1");378 gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B1"); 323 379 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE), 324 gr, l2, u, c, s1, true,5970, "#B2");380 gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B2"); 325 381 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH), 326 gr, l2, u, c, s1, true,5970, "#B3");382 gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B3"); 327 383 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST), 328 gr, l2, u, c, s1, true,5970, "#B4");384 gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B4"); 329 385 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST), 330 gr, l2, u, c, s1, true,5970, "#B5");386 gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B5"); 331 387 } 332 388 -
test/mip_test.cc
r674 r678 19 19 #include "test_tools.h" 20 20 21 #ifdef HAVE_CONFIG_H22 21 #include <lemon/config.h> 23 #endif24 22 25 23 #ifdef LEMON_HAVE_CPLEX -
tools/dimacs-solver.cc
r674 r691 120 120 ti.restart(); 121 121 NetworkSimplex<Digraph, Value> ns(g); 122 ns.lowerMap(lower). capacityMap(cap).costMap(cost).supplyMap(sup);123 if (sum_sup > 0) ns. problemType(ns.LEQ);122 ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup); 123 if (sum_sup > 0) ns.supplyType(ns.LEQ); 124 124 if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n'; 125 125 ti.restart();
Note: See TracChangeset
for help on using the changeset viewer.