# HG changeset patch
# User Alpar Juttner <alpar@cs.elte.hu>
# Date 1241029334 -3600
# Node ID e01957e96c6796a11d9f3bbf751eea8ae274221c
# Parent cb8270a9866086c20a37f916882b135f65c1eee6# Parent 8d289c89d43e63a309fc8a9a58100ecda310fb6f
Merge
diff -r cb8270a98660 -r e01957e96c67 CMakeLists.txt
--- a/CMakeLists.txt Wed Apr 29 16:15:29 2009 +0100
+++ b/CMakeLists.txt Wed Apr 29 19:22:14 2009 +0100
@@ -17,8 +17,6 @@
FIND_PACKAGE(CPLEX)
FIND_PACKAGE(COIN)
-ADD_DEFINITIONS(-DHAVE_CONFIG_H)
-
IF(MSVC)
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4250 /wd4355 /wd4800 /wd4996")
# Suppressed warnings:
@@ -28,8 +26,6 @@
# C4996: 'function': was declared deprecated
ENDIF(MSVC)
-ADD_DEFINITIONS(-DHAVE_CONFIG_H)
-
INCLUDE(CheckTypeSize)
CHECK_TYPE_SIZE("long long" LEMON_LONG_LONG)
diff -r cb8270a98660 -r e01957e96c67 Makefile.am
--- a/Makefile.am Wed Apr 29 16:15:29 2009 +0100
+++ b/Makefile.am Wed Apr 29 19:22:14 2009 +0100
@@ -11,11 +11,12 @@
m4/lx_check_cplex.m4 \
m4/lx_check_glpk.m4 \
m4/lx_check_soplex.m4 \
- m4/lx_check_clp.m4 \
- m4/lx_check_cbc.m4 \
+ m4/lx_check_coin.m4 \
CMakeLists.txt \
cmake/FindGhostscript.cmake \
+ cmake/FindCPLEX.cmake \
cmake/FindGLPK.cmake \
+ cmake/FindCOIN.cmake \
cmake/version.cmake.in \
cmake/version.cmake \
cmake/nsis/lemon.ico \
diff -r cb8270a98660 -r e01957e96c67 cmake/FindCOIN.cmake
--- a/cmake/FindCOIN.cmake Wed Apr 29 16:15:29 2009 +0100
+++ b/cmake/FindCOIN.cmake Wed Apr 29 19:22:14 2009 +0100
@@ -1,28 +1,48 @@
SET(COIN_ROOT_DIR "" CACHE PATH "COIN root directory")
FIND_PATH(COIN_INCLUDE_DIR coin/CoinUtilsConfig.h
- PATHS ${COIN_ROOT_DIR}/include)
-
-FIND_LIBRARY(COIN_CBC_LIBRARY libCbc
- PATHS ${COIN_ROOT_DIR}/lib)
-FIND_LIBRARY(COIN_CBC_SOLVER_LIBRARY libCbcSolver
- PATHS ${COIN_ROOT_DIR}/lib)
-FIND_LIBRARY(COIN_CGL_LIBRARY libCgl
- PATHS ${COIN_ROOT_DIR}/lib)
-FIND_LIBRARY(COIN_CLP_LIBRARY libClp
- PATHS ${COIN_ROOT_DIR}/lib)
-FIND_LIBRARY(COIN_COIN_UTILS_LIBRARY libCoinUtils
- PATHS ${COIN_ROOT_DIR}/lib)
-FIND_LIBRARY(COIN_OSI_LIBRARY libOsi
- PATHS ${COIN_ROOT_DIR}/lib)
-FIND_LIBRARY(COIN_OSI_CBC_LIBRARY libOsiCbc
- PATHS ${COIN_ROOT_DIR}/lib)
-FIND_LIBRARY(COIN_OSI_CLP_LIBRARY libOsiClp
- PATHS ${COIN_ROOT_DIR}/lib)
-FIND_LIBRARY(COIN_OSI_VOL_LIBRARY libOsiVol
- PATHS ${COIN_ROOT_DIR}/lib)
-FIND_LIBRARY(COIN_VOL_LIBRARY libVol
- PATHS ${COIN_ROOT_DIR}/lib)
+ HINTS ${COIN_ROOT_DIR}/include
+)
+FIND_LIBRARY(COIN_CBC_LIBRARY
+ NAMES Cbc libCbc
+ HINTS ${COIN_ROOT_DIR}/lib
+)
+FIND_LIBRARY(COIN_CBC_SOLVER_LIBRARY
+ NAMES CbcSolver libCbcSolver
+ HINTS ${COIN_ROOT_DIR}/lib
+)
+FIND_LIBRARY(COIN_CGL_LIBRARY
+ NAMES Cgl libCgl
+ HINTS ${COIN_ROOT_DIR}/lib
+)
+FIND_LIBRARY(COIN_CLP_LIBRARY
+ NAMES Clp libClp
+ HINTS ${COIN_ROOT_DIR}/lib
+)
+FIND_LIBRARY(COIN_COIN_UTILS_LIBRARY
+ NAMES CoinUtils libCoinUtils
+ HINTS ${COIN_ROOT_DIR}/lib
+)
+FIND_LIBRARY(COIN_OSI_LIBRARY
+ NAMES Osi libOsi
+ HINTS ${COIN_ROOT_DIR}/lib
+)
+FIND_LIBRARY(COIN_OSI_CBC_LIBRARY
+ NAMES OsiCbc libOsiCbc
+ HINTS ${COIN_ROOT_DIR}/lib
+)
+FIND_LIBRARY(COIN_OSI_CLP_LIBRARY
+ NAMES OsiClp libOsiClp
+ HINTS ${COIN_ROOT_DIR}/lib
+)
+FIND_LIBRARY(COIN_OSI_VOL_LIBRARY
+ NAMES OsiVol libOsiVol
+ HINTS ${COIN_ROOT_DIR}/lib
+)
+FIND_LIBRARY(COIN_VOL_LIBRARY
+ NAMES Vol libVol
+ HINTS ${COIN_ROOT_DIR}/lib
+)
INCLUDE(FindPackageHandleStandardArgs)
FIND_PACKAGE_HANDLE_STANDARD_ARGS(COIN DEFAULT_MSG
diff -r cb8270a98660 -r e01957e96c67 cmake/FindCPLEX.cmake
--- a/cmake/FindCPLEX.cmake Wed Apr 29 16:15:29 2009 +0100
+++ b/cmake/FindCPLEX.cmake Wed Apr 29 19:22:14 2009 +0100
@@ -1,21 +1,32 @@
+SET(CPLEX_ROOT_DIR "" CACHE PATH "CPLEX root directory")
+
FIND_PATH(CPLEX_INCLUDE_DIR
ilcplex/cplex.h
- PATHS "C:/ILOG/CPLEX91/include")
-
+ PATHS "C:/ILOG/CPLEX91/include"
+ PATHS "/opt/ilog/cplex91/include"
+ HINTS ${CPLEX_ROOT_DIR}/include
+)
FIND_LIBRARY(CPLEX_LIBRARY
- NAMES cplex91
- PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda")
+ cplex91
+ PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda"
+ PATHS "/opt/ilog/cplex91/bin"
+ HINTS ${CPLEX_ROOT_DIR}/bin
+)
INCLUDE(FindPackageHandleStandardArgs)
FIND_PACKAGE_HANDLE_STANDARD_ARGS(CPLEX DEFAULT_MSG CPLEX_LIBRARY CPLEX_INCLUDE_DIR)
FIND_PATH(CPLEX_BIN_DIR
cplex91.dll
- PATHS "C:/ILOG/CPLEX91/bin/x86_win32")
+ PATHS "C:/ILOG/CPLEX91/bin/x86_win32"
+)
IF(CPLEX_FOUND)
SET(CPLEX_INCLUDE_DIRS ${CPLEX_INCLUDE_DIR})
SET(CPLEX_LIBRARIES ${CPLEX_LIBRARY})
+ IF(CMAKE_SYSTEM_NAME STREQUAL "Linux")
+ SET(CPLEX_LIBRARIES "${CPLEX_LIBRARIES};m;pthread")
+ ENDIF(CMAKE_SYSTEM_NAME STREQUAL "Linux")
ENDIF(CPLEX_FOUND)
MARK_AS_ADVANCED(CPLEX_LIBRARY CPLEX_INCLUDE_DIR CPLEX_BIN_DIR)
diff -r cb8270a98660 -r e01957e96c67 cmake/FindGLPK.cmake
--- a/cmake/FindGLPK.cmake Wed Apr 29 16:15:29 2009 +0100
+++ b/cmake/FindGLPK.cmake Wed Apr 29 19:22:14 2009 +0100
@@ -1,16 +1,50 @@
+SET(GLPK_ROOT_DIR "" CACHE PATH "GLPK root directory")
+
SET(GLPK_REGKEY "[HKEY_LOCAL_MACHINE\\SOFTWARE\\GnuWin32\\Glpk;InstallPath]")
GET_FILENAME_COMPONENT(GLPK_ROOT_PATH ${GLPK_REGKEY} ABSOLUTE)
FIND_PATH(GLPK_INCLUDE_DIR
glpk.h
- PATHS ${GLPK_REGKEY}/include)
+ PATHS ${GLPK_REGKEY}/include
+ HINTS ${GLPK_ROOT_DIR}/include
+)
+FIND_LIBRARY(GLPK_LIBRARY
+ glpk
+ PATHS ${GLPK_REGKEY}/lib
+ HINTS ${GLPK_ROOT_DIR}/lib
+)
-FIND_LIBRARY(GLPK_LIBRARY
- NAMES glpk
- PATHS ${GLPK_REGKEY}/lib)
+IF(GLPK_INCLUDE_DIR AND GLPK_LIBRARY)
+ FILE(READ ${GLPK_INCLUDE_DIR}/glpk.h GLPK_GLPK_H)
+
+ STRING(REGEX MATCH "define[ ]+GLP_MAJOR_VERSION[ ]+[0-9]+" GLPK_MAJOR_VERSION_LINE "${GLPK_GLPK_H}")
+ STRING(REGEX REPLACE "define[ ]+GLP_MAJOR_VERSION[ ]+([0-9]+)" "\\1" GLPK_VERSION_MAJOR "${GLPK_MAJOR_VERSION_LINE}")
+
+ STRING(REGEX MATCH "define[ ]+GLP_MINOR_VERSION[ ]+[0-9]+" GLPK_MINOR_VERSION_LINE "${GLPK_GLPK_H}")
+ STRING(REGEX REPLACE "define[ ]+GLP_MINOR_VERSION[ ]+([0-9]+)" "\\1" GLPK_VERSION_MINOR "${GLPK_MINOR_VERSION_LINE}")
+
+ SET(GLPK_VERSION_STRING "${GLPK_VERSION_MAJOR}.${GLPK_VERSION_MINOR}")
+
+ IF(GLPK_FIND_VERSION)
+ IF(GLPK_FIND_VERSION_COUNT GREATER 2)
+ MESSAGE(SEND_ERROR "unexpected version string")
+ ENDIF(GLPK_FIND_VERSION_COUNT GREATER 2)
+
+ MATH(EXPR GLPK_REQUESTED_VERSION "${GLPK_FIND_VERSION_MAJOR}*100 + ${GLPK_FIND_VERSION_MINOR}")
+ MATH(EXPR GLPK_FOUND_VERSION "${GLPK_VERSION_MAJOR}*100 + ${GLPK_VERSION_MINOR}")
+
+ IF(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION)
+ SET(GLPK_PROPER_VERSION_FOUND FALSE)
+ ELSE(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION)
+ SET(GLPK_PROPER_VERSION_FOUND TRUE)
+ ENDIF(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION)
+ ELSE(GLPK_FIND_VERSION)
+ SET(GLPK_PROPER_VERSION_FOUND TRUE)
+ ENDIF(GLPK_FIND_VERSION)
+ENDIF(GLPK_INCLUDE_DIR AND GLPK_LIBRARY)
INCLUDE(FindPackageHandleStandardArgs)
-FIND_PACKAGE_HANDLE_STANDARD_ARGS(GLPK DEFAULT_MSG GLPK_LIBRARY GLPK_INCLUDE_DIR)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(GLPK DEFAULT_MSG GLPK_LIBRARY GLPK_INCLUDE_DIR GLPK_PROPER_VERSION_FOUND)
IF(GLPK_FOUND)
SET(GLPK_INCLUDE_DIRS ${GLPK_INCLUDE_DIR})
diff -r cb8270a98660 -r e01957e96c67 doc/groups.dox
--- a/doc/groups.dox Wed Apr 29 16:15:29 2009 +0100
+++ b/doc/groups.dox Wed Apr 29 19:22:14 2009 +0100
@@ -352,17 +352,17 @@
minimum total cost from a set of supply nodes to a set of demand nodes
in a network with capacity constraints (lower and upper bounds)
and arc costs.
-Formally, let \f$G=(V,A)\f$ be a digraph,
-\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
+Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
+\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
upper bounds for the flow values on the arcs, for which
-\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
-\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
-on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
+\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
+\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
+on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
signed supply values of the nodes.
If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
\f$-sup(u)\f$ demand.
-A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
+A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
of the following optimization problem.
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
@@ -404,7 +404,7 @@
The dual solution of the minimum cost flow problem is represented by node
potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
-An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
+An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
node potentials the following \e complementary \e slackness optimality
conditions hold.
@@ -413,15 +413,15 @@
- if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
- if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
- if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
- - For all \f$u\in V\f$:
+ - For all \f$u\in V\f$ nodes:
- if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
then \f$\pi(u)=0\f$.
Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
-\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
+\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
-All algorithms provide dual solution (node potentials) as well
+All algorithms provide dual solution (node potentials) as well,
if an optimal flow is found.
LEMON contains several algorithms for solving minimum cost flow problems.
diff -r cb8270a98660 -r e01957e96c67 lemon/Makefile.am
--- a/lemon/Makefile.am Wed Apr 29 16:15:29 2009 +0100
+++ b/lemon/Makefile.am Wed Apr 29 19:22:14 2009 +0100
@@ -15,7 +15,8 @@
lemon/random.cc \
lemon/bits/windows.cc
-
+nodist_lemon_HEADERS = lemon/config.h
+
lemon_libemon_la_CXXFLAGS = \
$(AM_CXXFLAGS) \
$(GLPK_CFLAGS) \
@@ -57,11 +58,11 @@
lemon/assert.h \
lemon/bfs.h \
lemon/bin_heap.h \
+ lemon/cbc.h \
lemon/circulation.h \
lemon/clp.h \
lemon/color.h \
lemon/concept_check.h \
- lemon/config.h \
lemon/connectivity.h \
lemon/counter.h \
lemon/core.h \
diff -r cb8270a98660 -r e01957e96c67 lemon/circulation.h
--- a/lemon/circulation.h Wed Apr 29 16:15:29 2009 +0100
+++ b/lemon/circulation.h Wed Apr 29 19:22:14 2009 +0100
@@ -64,15 +64,15 @@
/// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
typedef SM SupplyMap;
- /// \brief The type of the flow values.
- typedef typename SupplyMap::Value Flow;
+ /// \brief The type of the flow and supply values.
+ typedef typename SupplyMap::Value Value;
/// \brief The type of the map that stores the flow values.
///
/// The type of the map that stores the flow values.
/// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
/// concept.
- typedef typename Digraph::template ArcMap<Flow> FlowMap;
+ typedef typename Digraph::template ArcMap<Value> FlowMap;
/// \brief Instantiates a FlowMap.
///
@@ -104,7 +104,7 @@
/// \brief The tolerance used by the algorithm
///
/// The tolerance used by the algorithm to handle inexact computation.
- typedef lemon::Tolerance<Flow> Tolerance;
+ typedef lemon::Tolerance<Value> Tolerance;
};
@@ -187,8 +187,8 @@
typedef TR Traits;
///The type of the digraph the algorithm runs on.
typedef typename Traits::Digraph Digraph;
- ///The type of the flow values.
- typedef typename Traits::Flow Flow;
+ ///The type of the flow and supply values.
+ typedef typename Traits::Value Value;
///The type of the lower bound map.
typedef typename Traits::LowerMap LowerMap;
@@ -221,7 +221,7 @@
Elevator* _level;
bool _local_level;
- typedef typename Digraph::template NodeMap<Flow> ExcessMap;
+ typedef typename Digraph::template NodeMap<Value> ExcessMap;
ExcessMap* _excess;
Tolerance _tol;
@@ -530,7 +530,7 @@
(*_excess)[_g.target(e)] += (*_lo)[e];
(*_excess)[_g.source(e)] -= (*_lo)[e];
} else {
- Flow fc = -(*_excess)[_g.target(e)];
+ Value fc = -(*_excess)[_g.target(e)];
_flow->set(e, fc);
(*_excess)[_g.target(e)] = 0;
(*_excess)[_g.source(e)] -= fc;
@@ -563,11 +563,11 @@
while((act=_level->highestActive())!=INVALID) {
int actlevel=(*_level)[act];
int mlevel=_node_num;
- Flow exc=(*_excess)[act];
+ Value exc=(*_excess)[act];
for(OutArcIt e(_g,act);e!=INVALID; ++e) {
Node v = _g.target(e);
- Flow fc=(*_up)[e]-(*_flow)[e];
+ Value fc=(*_up)[e]-(*_flow)[e];
if(!_tol.positive(fc)) continue;
if((*_level)[v]<actlevel) {
if(!_tol.less(fc, exc)) {
@@ -591,7 +591,7 @@
}
for(InArcIt e(_g,act);e!=INVALID; ++e) {
Node v = _g.source(e);
- Flow fc=(*_flow)[e]-(*_lo)[e];
+ Value fc=(*_flow)[e]-(*_lo)[e];
if(!_tol.positive(fc)) continue;
if((*_level)[v]<actlevel) {
if(!_tol.less(fc, exc)) {
@@ -661,13 +661,13 @@
///@{
- /// \brief Returns the flow on the given arc.
+ /// \brief Returns the flow value on the given arc.
///
- /// Returns the flow on the given arc.
+ /// Returns the flow value on the given arc.
///
/// \pre Either \ref run() or \ref init() must be called before
/// using this function.
- Flow flow(const Arc& arc) const {
+ Value flow(const Arc& arc) const {
return (*_flow)[arc];
}
@@ -750,7 +750,7 @@
if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
for(NodeIt n(_g);n!=INVALID;++n)
{
- Flow dif=-(*_supply)[n];
+ Value dif=-(*_supply)[n];
for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
if(_tol.negative(dif)) return false;
@@ -765,10 +765,10 @@
///\sa barrierMap()
bool checkBarrier() const
{
- Flow delta=0;
- Flow inf_cap = std::numeric_limits<Flow>::has_infinity ?
- std::numeric_limits<Flow>::infinity() :
- std::numeric_limits<Flow>::max();
+ Value delta=0;
+ Value inf_cap = std::numeric_limits<Value>::has_infinity ?
+ std::numeric_limits<Value>::infinity() :
+ std::numeric_limits<Value>::max();
for(NodeIt n(_g);n!=INVALID;++n)
if(barrier(n))
delta-=(*_supply)[n];
diff -r cb8270a98660 -r e01957e96c67 lemon/core.h
--- a/lemon/core.h Wed Apr 29 16:15:29 2009 +0100
+++ b/lemon/core.h Wed Apr 29 19:22:14 2009 +0100
@@ -22,7 +22,7 @@
#include <vector>
#include <algorithm>
-#include <lemon/core.h>
+#include <lemon/config.h>
#include <lemon/bits/enable_if.h>
#include <lemon/bits/traits.h>
#include <lemon/assert.h>
diff -r cb8270a98660 -r e01957e96c67 lemon/network_simplex.h
--- a/lemon/network_simplex.h Wed Apr 29 16:15:29 2009 +0100
+++ b/lemon/network_simplex.h Wed Apr 29 19:22:14 2009 +0100
@@ -30,9 +30,6 @@
#include <lemon/core.h>
#include <lemon/math.h>
-#include <lemon/maps.h>
-#include <lemon/circulation.h>
-#include <lemon/adaptors.h>
namespace lemon {
@@ -50,14 +47,19 @@
///
/// In general this class is the fastest implementation available
/// in LEMON for the minimum cost flow problem.
- /// Moreover it supports both direction of the supply/demand inequality
- /// constraints. For more information see \ref ProblemType.
+ /// Moreover it supports both directions of the supply/demand inequality
+ /// constraints. For more information see \ref SupplyType.
+ ///
+ /// Most of the parameters of the problem (except for the digraph)
+ /// can be given using separate functions, and the algorithm can be
+ /// executed using the \ref run() function. If some parameters are not
+ /// specified, then default values will be used.
///
/// \tparam GR The digraph type the algorithm runs on.
- /// \tparam F The value type used for flow amounts, capacity bounds
+ /// \tparam V The value type used for flow amounts, capacity bounds
/// and supply values in the algorithm. By default it is \c int.
/// \tparam C The value type used for costs and potentials in the
- /// algorithm. By default it is the same as \c F.
+ /// algorithm. By default it is the same as \c V.
///
/// \warning Both value types must be signed and all input data must
/// be integer.
@@ -65,34 +67,92 @@
/// \note %NetworkSimplex provides five different pivot rule
/// implementations, from which the most efficient one is used
/// by default. For more information see \ref PivotRule.
- template <typename GR, typename F = int, typename C = F>
+ template <typename GR, typename V = int, typename C = V>
class NetworkSimplex
{
public:
- /// The flow type of the algorithm
- typedef F Flow;
- /// The cost type of the algorithm
+ /// The type of the flow amounts, capacity bounds and supply values
+ typedef V Value;
+ /// The type of the arc costs
typedef C Cost;
-#ifdef DOXYGEN
- /// The type of the flow map
- typedef GR::ArcMap<Flow> FlowMap;
- /// The type of the potential map
- typedef GR::NodeMap<Cost> PotentialMap;
-#else
- /// The type of the flow map
- typedef typename GR::template ArcMap<Flow> FlowMap;
- /// The type of the potential map
- typedef typename GR::template NodeMap<Cost> PotentialMap;
-#endif
public:
- /// \brief Enum type for selecting the pivot rule.
+ /// \brief Problem type constants for the \c run() function.
///
- /// Enum type for selecting the pivot rule for the \ref run()
+ /// Enum type containing the problem type constants that can be
+ /// returned by the \ref run() function of the algorithm.
+ enum ProblemType {
+ /// The problem has no feasible solution (flow).
+ INFEASIBLE,
+ /// The problem has optimal solution (i.e. it is feasible and
+ /// bounded), and the algorithm has found optimal flow and node
+ /// potentials (primal and dual solutions).
+ OPTIMAL,
+ /// The objective function of the problem is unbounded, i.e.
+ /// there is a directed cycle having negative total cost and
+ /// infinite upper bound.
+ UNBOUNDED
+ };
+
+ /// \brief Constants for selecting the type of the supply constraints.
+ ///
+ /// Enum type containing constants for selecting the supply type,
+ /// i.e. the direction of the inequalities in the supply/demand
+ /// constraints of the \ref min_cost_flow "minimum cost flow problem".
+ ///
+ /// The default supply type is \c GEQ, since this form is supported
+ /// by other minimum cost flow algorithms and the \ref Circulation
+ /// algorithm, as well.
+ /// The \c LEQ problem type can be selected using the \ref supplyType()
/// function.
///
+ /// Note that the equality form is a special case of both supply types.
+ enum SupplyType {
+
+ /// This option means that there are <em>"greater or equal"</em>
+ /// supply/demand constraints in the definition, i.e. the exact
+ /// formulation of the problem is the following.
+ /**
+ \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
+ \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
+ sup(u) \quad \forall u\in V \f]
+ \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
+ */
+ /// It means that the total demand must be greater or equal to the
+ /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
+ /// negative) and all the supplies have to be carried out from
+ /// the supply nodes, but there could be demands that are not
+ /// satisfied.
+ GEQ,
+ /// It is just an alias for the \c GEQ option.
+ CARRY_SUPPLIES = GEQ,
+
+ /// This option means that there are <em>"less or equal"</em>
+ /// supply/demand constraints in the definition, i.e. the exact
+ /// formulation of the problem is the following.
+ /**
+ \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
+ \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
+ sup(u) \quad \forall u\in V \f]
+ \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
+ */
+ /// It means that the total demand must be less or equal to the
+ /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
+ /// positive) and all the demands have to be satisfied, but there
+ /// could be supplies that are not carried out from the supply
+ /// nodes.
+ LEQ,
+ /// It is just an alias for the \c LEQ option.
+ SATISFY_DEMANDS = LEQ
+ };
+
+ /// \brief Constants for selecting the pivot rule.
+ ///
+ /// Enum type containing constants for selecting the pivot rule for
+ /// the \ref run() function.
+ ///
/// \ref NetworkSimplex provides five different pivot rule
/// implementations that significantly affect the running time
/// of the algorithm.
@@ -131,71 +191,15 @@
ALTERING_LIST
};
- /// \brief Enum type for selecting the problem type.
- ///
- /// Enum type for selecting the problem type, i.e. the direction of
- /// the inequalities in the supply/demand constraints of the
- /// \ref min_cost_flow "minimum cost flow problem".
- ///
- /// The default problem type is \c GEQ, since this form is supported
- /// by other minimum cost flow algorithms and the \ref Circulation
- /// algorithm as well.
- /// The \c LEQ problem type can be selected using the \ref problemType()
- /// function.
- ///
- /// Note that the equality form is a special case of both problem type.
- enum ProblemType {
-
- /// This option means that there are "<em>greater or equal</em>"
- /// constraints in the defintion, i.e. the exact formulation of the
- /// problem is the following.
- /**
- \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
- \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
- sup(u) \quad \forall u\in V \f]
- \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
- */
- /// It means that the total demand must be greater or equal to the
- /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
- /// negative) and all the supplies have to be carried out from
- /// the supply nodes, but there could be demands that are not
- /// satisfied.
- GEQ,
- /// It is just an alias for the \c GEQ option.
- CARRY_SUPPLIES = GEQ,
-
- /// This option means that there are "<em>less or equal</em>"
- /// constraints in the defintion, i.e. the exact formulation of the
- /// problem is the following.
- /**
- \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
- \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
- sup(u) \quad \forall u\in V \f]
- \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
- */
- /// It means that the total demand must be less or equal to the
- /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
- /// positive) and all the demands have to be satisfied, but there
- /// could be supplies that are not carried out from the supply
- /// nodes.
- LEQ,
- /// It is just an alias for the \c LEQ option.
- SATISFY_DEMANDS = LEQ
- };
-
private:
TEMPLATE_DIGRAPH_TYPEDEFS(GR);
- typedef typename GR::template ArcMap<Flow> FlowArcMap;
- typedef typename GR::template ArcMap<Cost> CostArcMap;
- typedef typename GR::template NodeMap<Flow> FlowNodeMap;
-
typedef std::vector<Arc> ArcVector;
typedef std::vector<Node> NodeVector;
typedef std::vector<int> IntVector;
typedef std::vector<bool> BoolVector;
- typedef std::vector<Flow> FlowVector;
+ typedef std::vector<Value> ValueVector;
typedef std::vector<Cost> CostVector;
// State constants for arcs
@@ -213,32 +217,23 @@
int _arc_num;
// Parameters of the problem
- FlowArcMap *_plower;
- FlowArcMap *_pupper;
- CostArcMap *_pcost;
- FlowNodeMap *_psupply;
- bool _pstsup;
- Node _psource, _ptarget;
- Flow _pstflow;
- ProblemType _ptype;
-
- // Result maps
- FlowMap *_flow_map;
- PotentialMap *_potential_map;
- bool _local_flow;
- bool _local_potential;
+ bool _have_lower;
+ SupplyType _stype;
+ Value _sum_supply;
// Data structures for storing the digraph
IntNodeMap _node_id;
- ArcVector _arc_ref;
+ IntArcMap _arc_id;
IntVector _source;
IntVector _target;
// Node and arc data
- FlowVector _cap;
+ ValueVector _lower;
+ ValueVector _upper;
+ ValueVector _cap;
CostVector _cost;
- FlowVector _supply;
- FlowVector _flow;
+ ValueVector _supply;
+ ValueVector _flow;
CostVector _pi;
// Data for storing the spanning tree structure
@@ -257,7 +252,16 @@
int in_arc, join, u_in, v_in, u_out, v_out;
int first, second, right, last;
int stem, par_stem, new_stem;
- Flow delta;
+ Value delta;
+
+ public:
+
+ /// \brief Constant for infinite upper bounds (capacities).
+ ///
+ /// Constant for infinite upper bounds (capacities).
+ /// It is \c std::numeric_limits<Value>::infinity() if available,
+ /// \c std::numeric_limits<Value>::max() otherwise.
+ const Value INF;
private:
@@ -659,25 +663,68 @@
///
/// \param graph The digraph the algorithm runs on.
NetworkSimplex(const GR& graph) :
- _graph(graph),
- _plower(NULL), _pupper(NULL), _pcost(NULL),
- _psupply(NULL), _pstsup(false), _ptype(GEQ),
- _flow_map(NULL), _potential_map(NULL),
- _local_flow(false), _local_potential(false),
- _node_id(graph)
+ _graph(graph), _node_id(graph), _arc_id(graph),
+ INF(std::numeric_limits<Value>::has_infinity ?
+ std::numeric_limits<Value>::infinity() :
+ std::numeric_limits<Value>::max())
{
- LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
- std::numeric_limits<Flow>::is_signed,
- "The flow type of NetworkSimplex must be signed integer");
- LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
- std::numeric_limits<Cost>::is_signed,
- "The cost type of NetworkSimplex must be signed integer");
- }
+ // Check the value types
+ LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
+ "The flow type of NetworkSimplex must be signed");
+ LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
+ "The cost type of NetworkSimplex must be signed");
+
+ // Resize vectors
+ _node_num = countNodes(_graph);
+ _arc_num = countArcs(_graph);
+ int all_node_num = _node_num + 1;
+ int all_arc_num = _arc_num + _node_num;
- /// Destructor.
- ~NetworkSimplex() {
- if (_local_flow) delete _flow_map;
- if (_local_potential) delete _potential_map;
+ _source.resize(all_arc_num);
+ _target.resize(all_arc_num);
+
+ _lower.resize(all_arc_num);
+ _upper.resize(all_arc_num);
+ _cap.resize(all_arc_num);
+ _cost.resize(all_arc_num);
+ _supply.resize(all_node_num);
+ _flow.resize(all_arc_num);
+ _pi.resize(all_node_num);
+
+ _parent.resize(all_node_num);
+ _pred.resize(all_node_num);
+ _forward.resize(all_node_num);
+ _thread.resize(all_node_num);
+ _rev_thread.resize(all_node_num);
+ _succ_num.resize(all_node_num);
+ _last_succ.resize(all_node_num);
+ _state.resize(all_arc_num);
+
+ // Copy the graph (store the arcs in a mixed order)
+ int i = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+ _node_id[n] = i;
+ }
+ int k = std::max(int(std::sqrt(double(_arc_num))), 10);
+ i = 0;
+ for (ArcIt a(_graph); a != INVALID; ++a) {
+ _arc_id[a] = i;
+ _source[i] = _node_id[_graph.source(a)];
+ _target[i] = _node_id[_graph.target(a)];
+ if ((i += k) >= _arc_num) i = (i % k) + 1;
+ }
+
+ // Initialize maps
+ for (int i = 0; i != _node_num; ++i) {
+ _supply[i] = 0;
+ }
+ for (int i = 0; i != _arc_num; ++i) {
+ _lower[i] = 0;
+ _upper[i] = INF;
+ _cost[i] = 1;
+ }
+ _have_lower = false;
+ _stype = GEQ;
}
/// \name Parameters
@@ -689,21 +736,19 @@
/// \brief Set the lower bounds on the arcs.
///
/// This function sets the lower bounds on the arcs.
- /// If neither this function nor \ref boundMaps() is used before
- /// calling \ref run(), the lower bounds will be set to zero
- /// on all arcs.
+ /// If it is not used before calling \ref run(), the lower bounds
+ /// will be set to zero on all arcs.
///
/// \param map An arc map storing the lower bounds.
- /// Its \c Value type must be convertible to the \c Flow type
+ /// Its \c Value type must be convertible to the \c Value type
/// of the algorithm.
///
/// \return <tt>(*this)</tt>
- template <typename LOWER>
- NetworkSimplex& lowerMap(const LOWER& map) {
- delete _plower;
- _plower = new FlowArcMap(_graph);
+ template <typename LowerMap>
+ NetworkSimplex& lowerMap(const LowerMap& map) {
+ _have_lower = true;
for (ArcIt a(_graph); a != INVALID; ++a) {
- (*_plower)[a] = map[a];
+ _lower[_arc_id[a]] = map[a];
}
return *this;
}
@@ -711,63 +756,23 @@
/// \brief Set the upper bounds (capacities) on the arcs.
///
/// This function sets the upper bounds (capacities) on the arcs.
- /// If none of the functions \ref upperMap(), \ref capacityMap()
- /// and \ref boundMaps() is used before calling \ref run(),
- /// the upper bounds (capacities) will be set to
- /// \c std::numeric_limits<Flow>::max() on all arcs.
+ /// If it is not used before calling \ref run(), the upper bounds
+ /// will be set to \ref INF on all arcs (i.e. the flow value will be
+ /// unbounded from above on each arc).
///
/// \param map An arc map storing the upper bounds.
- /// Its \c Value type must be convertible to the \c Flow type
+ /// Its \c Value type must be convertible to the \c Value type
/// of the algorithm.
///
/// \return <tt>(*this)</tt>
- template<typename UPPER>
- NetworkSimplex& upperMap(const UPPER& map) {
- delete _pupper;
- _pupper = new FlowArcMap(_graph);
+ template<typename UpperMap>
+ NetworkSimplex& upperMap(const UpperMap& map) {
for (ArcIt a(_graph); a != INVALID; ++a) {
- (*_pupper)[a] = map[a];
+ _upper[_arc_id[a]] = map[a];
}
return *this;
}
- /// \brief Set the upper bounds (capacities) on the arcs.
- ///
- /// This function sets the upper bounds (capacities) on the arcs.
- /// It is just an alias for \ref upperMap().
- ///
- /// \return <tt>(*this)</tt>
- template<typename CAP>
- NetworkSimplex& capacityMap(const CAP& map) {
- return upperMap(map);
- }
-
- /// \brief Set the lower and upper bounds on the arcs.
- ///
- /// This function sets the lower and upper bounds on the arcs.
- /// If neither this function nor \ref lowerMap() is used before
- /// calling \ref run(), the lower bounds will be set to zero
- /// on all arcs.
- /// If none of the functions \ref upperMap(), \ref capacityMap()
- /// and \ref boundMaps() is used before calling \ref run(),
- /// the upper bounds (capacities) will be set to
- /// \c std::numeric_limits<Flow>::max() on all arcs.
- ///
- /// \param lower An arc map storing the lower bounds.
- /// \param upper An arc map storing the upper bounds.
- ///
- /// The \c Value type of the maps must be convertible to the
- /// \c Flow type of the algorithm.
- ///
- /// \note This function is just a shortcut of calling \ref lowerMap()
- /// and \ref upperMap() separately.
- ///
- /// \return <tt>(*this)</tt>
- template <typename LOWER, typename UPPER>
- NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
- return lowerMap(lower).upperMap(upper);
- }
-
/// \brief Set the costs of the arcs.
///
/// This function sets the costs of the arcs.
@@ -779,12 +784,10 @@
/// of the algorithm.
///
/// \return <tt>(*this)</tt>
- template<typename COST>
- NetworkSimplex& costMap(const COST& map) {
- delete _pcost;
- _pcost = new CostArcMap(_graph);
+ template<typename CostMap>
+ NetworkSimplex& costMap(const CostMap& map) {
for (ArcIt a(_graph); a != INVALID; ++a) {
- (*_pcost)[a] = map[a];
+ _cost[_arc_id[a]] = map[a];
}
return *this;
}
@@ -797,17 +800,14 @@
/// (It makes sense only if non-zero lower bounds are given.)
///
/// \param map A node map storing the supply values.
- /// Its \c Value type must be convertible to the \c Flow type
+ /// Its \c Value type must be convertible to the \c Value type
/// of the algorithm.
///
/// \return <tt>(*this)</tt>
- template<typename SUP>
- NetworkSimplex& supplyMap(const SUP& map) {
- delete _psupply;
- _pstsup = false;
- _psupply = new FlowNodeMap(_graph);
+ template<typename SupplyMap>
+ NetworkSimplex& supplyMap(const SupplyMap& map) {
for (NodeIt n(_graph); n != INVALID; ++n) {
- (*_psupply)[n] = map[n];
+ _supply[_node_id[n]] = map[n];
}
return *this;
}
@@ -820,71 +820,39 @@
/// calling \ref run(), the supply of each node will be set to zero.
/// (It makes sense only if non-zero lower bounds are given.)
///
+ /// Using this function has the same effect as using \ref supplyMap()
+ /// with such a map in which \c k is assigned to \c s, \c -k is
+ /// assigned to \c t and all other nodes have zero supply value.
+ ///
/// \param s The source node.
/// \param t The target node.
/// \param k The required amount of flow from node \c s to node \c t
/// (i.e. the supply of \c s and the demand of \c t).
///
/// \return <tt>(*this)</tt>
- NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
- delete _psupply;
- _psupply = NULL;
- _pstsup = true;
- _psource = s;
- _ptarget = t;
- _pstflow = k;
+ NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
+ for (int i = 0; i != _node_num; ++i) {
+ _supply[i] = 0;
+ }
+ _supply[_node_id[s]] = k;
+ _supply[_node_id[t]] = -k;
return *this;
}
- /// \brief Set the problem type.
+ /// \brief Set the type of the supply constraints.
///
- /// This function sets the problem type for the algorithm.
- /// If it is not used before calling \ref run(), the \ref GEQ problem
+ /// This function sets the type of the supply/demand constraints.
+ /// If it is not used before calling \ref run(), the \ref GEQ supply
/// type will be used.
///
- /// For more information see \ref ProblemType.
+ /// For more information see \ref SupplyType.
///
/// \return <tt>(*this)</tt>
- NetworkSimplex& problemType(ProblemType problem_type) {
- _ptype = problem_type;
+ NetworkSimplex& supplyType(SupplyType supply_type) {
+ _stype = supply_type;
return *this;
}
- /// \brief Set the flow map.
- ///
- /// This function sets the flow map.
- /// If it is not used before calling \ref run(), an instance will
- /// be allocated automatically. The destructor deallocates this
- /// automatically allocated map, of course.
- ///
- /// \return <tt>(*this)</tt>
- NetworkSimplex& flowMap(FlowMap& map) {
- if (_local_flow) {
- delete _flow_map;
- _local_flow = false;
- }
- _flow_map = ↦
- return *this;
- }
-
- /// \brief Set the potential map.
- ///
- /// This function sets the potential map, which is used for storing
- /// the dual solution.
- /// If it is not used before calling \ref run(), an instance will
- /// be allocated automatically. The destructor deallocates this
- /// automatically allocated map, of course.
- ///
- /// \return <tt>(*this)</tt>
- NetworkSimplex& potentialMap(PotentialMap& map) {
- if (_local_potential) {
- delete _potential_map;
- _local_potential = false;
- }
- _potential_map = ↦
- return *this;
- }
-
/// @}
/// \name Execution Control
@@ -896,13 +864,12 @@
///
/// This function runs the algorithm.
/// The paramters can be specified using functions \ref lowerMap(),
- /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
- /// \ref costMap(), \ref supplyMap(), \ref stSupply(),
- /// \ref problemType(), \ref flowMap() and \ref potentialMap().
+ /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
+ /// \ref supplyType().
/// For example,
/// \code
/// NetworkSimplex<ListDigraph> ns(graph);
- /// ns.boundMaps(lower, upper).costMap(cost)
+ /// ns.lowerMap(lower).upperMap(upper).costMap(cost)
/// .supplyMap(sup).run();
/// \endcode
///
@@ -910,33 +877,44 @@
/// that have been given are kept for the next call, unless
/// \ref reset() is called, thus only the modified parameters
/// have to be set again. See \ref reset() for examples.
+ /// However the underlying digraph must not be modified after this
+ /// class have been constructed, since it copies and extends the graph.
///
/// \param pivot_rule The pivot rule that will be used during the
/// algorithm. For more information see \ref PivotRule.
///
- /// \return \c true if a feasible flow can be found.
- bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
- return init() && start(pivot_rule);
+ /// \return \c INFEASIBLE if no feasible flow exists,
+ /// \n \c OPTIMAL if the problem has optimal solution
+ /// (i.e. it is feasible and bounded), and the algorithm has found
+ /// optimal flow and node potentials (primal and dual solutions),
+ /// \n \c UNBOUNDED if the objective function of the problem is
+ /// unbounded, i.e. there is a directed cycle having negative total
+ /// cost and infinite upper bound.
+ ///
+ /// \see ProblemType, PivotRule
+ ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
+ if (!init()) return INFEASIBLE;
+ return start(pivot_rule);
}
/// \brief Reset all the parameters that have been given before.
///
/// This function resets all the paramaters that have been given
/// before using functions \ref lowerMap(), \ref upperMap(),
- /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
- /// \ref supplyMap(), \ref stSupply(), \ref problemType(),
- /// \ref flowMap() and \ref potentialMap().
+ /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
///
/// It is useful for multiple run() calls. If this function is not
/// used, all the parameters given before are kept for the next
/// \ref run() call.
+ /// However the underlying digraph must not be modified after this
+ /// class have been constructed, since it copies and extends the graph.
///
/// For example,
/// \code
/// NetworkSimplex<ListDigraph> ns(graph);
///
/// // First run
- /// ns.lowerMap(lower).capacityMap(cap).costMap(cost)
+ /// ns.lowerMap(lower).upperMap(upper).costMap(cost)
/// .supplyMap(sup).run();
///
/// // Run again with modified cost map (reset() is not called,
@@ -947,29 +925,22 @@
/// // Run again from scratch using reset()
/// // (the lower bounds will be set to zero on all arcs)
/// ns.reset();
- /// ns.capacityMap(cap).costMap(cost)
+ /// ns.upperMap(capacity).costMap(cost)
/// .supplyMap(sup).run();
/// \endcode
///
/// \return <tt>(*this)</tt>
NetworkSimplex& reset() {
- delete _plower;
- delete _pupper;
- delete _pcost;
- delete _psupply;
- _plower = NULL;
- _pupper = NULL;
- _pcost = NULL;
- _psupply = NULL;
- _pstsup = false;
- _ptype = GEQ;
- if (_local_flow) delete _flow_map;
- if (_local_potential) delete _potential_map;
- _flow_map = NULL;
- _potential_map = NULL;
- _local_flow = false;
- _local_potential = false;
-
+ for (int i = 0; i != _node_num; ++i) {
+ _supply[i] = 0;
+ }
+ for (int i = 0; i != _arc_num; ++i) {
+ _lower[i] = 0;
+ _upper[i] = INF;
+ _cost[i] = 1;
+ }
+ _have_lower = false;
+ _stype = GEQ;
return *this;
}
@@ -985,7 +956,7 @@
/// \brief Return the total cost of the found flow.
///
/// This function returns the total cost of the found flow.
- /// The complexity of the function is O(e).
+ /// Its complexity is O(e).
///
/// \note The return type of the function can be specified as a
/// template parameter. For example,
@@ -997,15 +968,12 @@
/// function.
///
/// \pre \ref run() must be called before using this function.
- template <typename Num>
- Num totalCost() const {
- Num c = 0;
- if (_pcost) {
- for (ArcIt e(_graph); e != INVALID; ++e)
- c += (*_flow_map)[e] * (*_pcost)[e];
- } else {
- for (ArcIt e(_graph); e != INVALID; ++e)
- c += (*_flow_map)[e];
+ template <typename Number>
+ Number totalCost() const {
+ Number c = 0;
+ for (ArcIt a(_graph); a != INVALID; ++a) {
+ int i = _arc_id[a];
+ c += Number(_flow[i]) * Number(_cost[i]);
}
return c;
}
@@ -1021,18 +989,22 @@
/// This function returns the flow on the given arc.
///
/// \pre \ref run() must be called before using this function.
- Flow flow(const Arc& a) const {
- return (*_flow_map)[a];
+ Value flow(const Arc& a) const {
+ return _flow[_arc_id[a]];
}
- /// \brief Return a const reference to the flow map.
+ /// \brief Return the flow map (the primal solution).
///
- /// This function returns a const reference to an arc map storing
- /// the found flow.
+ /// This function copies the flow value on each arc into the given
+ /// map. The \c Value type of the algorithm must be convertible to
+ /// the \c Value type of the map.
///
/// \pre \ref run() must be called before using this function.
- const FlowMap& flowMap() const {
- return *_flow_map;
+ template <typename FlowMap>
+ void flowMap(FlowMap &map) const {
+ for (ArcIt a(_graph); a != INVALID; ++a) {
+ map.set(a, _flow[_arc_id[a]]);
+ }
}
/// \brief Return the potential (dual value) of the given node.
@@ -1042,19 +1014,22 @@
///
/// \pre \ref run() must be called before using this function.
Cost potential(const Node& n) const {
- return (*_potential_map)[n];
+ return _pi[_node_id[n]];
}
- /// \brief Return a const reference to the potential map
- /// (the dual solution).
+ /// \brief Return the potential map (the dual solution).
///
- /// This function returns a const reference to a node map storing
- /// the found potentials, which form the dual solution of the
- /// \ref min_cost_flow "minimum cost flow" problem.
+ /// This function copies the potential (dual value) of each node
+ /// into the given map.
+ /// The \c Cost type of the algorithm must be convertible to the
+ /// \c Value type of the map.
///
/// \pre \ref run() must be called before using this function.
- const PotentialMap& potentialMap() const {
- return *_potential_map;
+ template <typename PotentialMap>
+ void potentialMap(PotentialMap &map) const {
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ map.set(n, _pi[_node_id[n]]);
+ }
}
/// @}
@@ -1063,245 +1038,82 @@
// Initialize internal data structures
bool init() {
- // Initialize result maps
- if (!_flow_map) {
- _flow_map = new FlowMap(_graph);
- _local_flow = true;
+ if (_node_num == 0) return false;
+
+ // Check the sum of supply values
+ _sum_supply = 0;
+ for (int i = 0; i != _node_num; ++i) {
+ _sum_supply += _supply[i];
}
- if (!_potential_map) {
- _potential_map = new PotentialMap(_graph);
- _local_potential = true;
+ if ( !((_stype == GEQ && _sum_supply <= 0) ||
+ (_stype == LEQ && _sum_supply >= 0)) ) return false;
+
+ // Remove non-zero lower bounds
+ if (_have_lower) {
+ for (int i = 0; i != _arc_num; ++i) {
+ Value c = _lower[i];
+ if (c >= 0) {
+ _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
+ } else {
+ _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
+ }
+ _supply[_source[i]] -= c;
+ _supply[_target[i]] += c;
+ }
+ } else {
+ for (int i = 0; i != _arc_num; ++i) {
+ _cap[i] = _upper[i];
+ }
}
- // Initialize vectors
- _node_num = countNodes(_graph);
- _arc_num = countArcs(_graph);
- int all_node_num = _node_num + 1;
- int all_arc_num = _arc_num + _node_num;
- if (_node_num == 0) return false;
-
- _arc_ref.resize(_arc_num);
- _source.resize(all_arc_num);
- _target.resize(all_arc_num);
-
- _cap.resize(all_arc_num);
- _cost.resize(all_arc_num);
- _supply.resize(all_node_num);
- _flow.resize(all_arc_num);
- _pi.resize(all_node_num);
-
- _parent.resize(all_node_num);
- _pred.resize(all_node_num);
- _forward.resize(all_node_num);
- _thread.resize(all_node_num);
- _rev_thread.resize(all_node_num);
- _succ_num.resize(all_node_num);
- _last_succ.resize(all_node_num);
- _state.resize(all_arc_num);
-
- // Initialize node related data
- bool valid_supply = true;
- Flow sum_supply = 0;
- if (!_pstsup && !_psupply) {
- _pstsup = true;
- _psource = _ptarget = NodeIt(_graph);
- _pstflow = 0;
- }
- if (_psupply) {
- int i = 0;
- for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
- _node_id[n] = i;
- _supply[i] = (*_psupply)[n];
- sum_supply += _supply[i];
+ // Initialize artifical cost
+ Cost ART_COST;
+ if (std::numeric_limits<Cost>::is_exact) {
+ ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
+ } else {
+ ART_COST = std::numeric_limits<Cost>::min();
+ for (int i = 0; i != _arc_num; ++i) {
+ if (_cost[i] > ART_COST) ART_COST = _cost[i];
}
- valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
- (_ptype == LEQ && sum_supply >= 0);
- } else {
- int i = 0;
- for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
- _node_id[n] = i;
- _supply[i] = 0;
- }
- _supply[_node_id[_psource]] = _pstflow;
- _supply[_node_id[_ptarget]] = -_pstflow;
- }
- if (!valid_supply) return false;
-
- // Infinite capacity value
- Flow inf_cap =
- std::numeric_limits<Flow>::has_infinity ?
- std::numeric_limits<Flow>::infinity() :
- std::numeric_limits<Flow>::max();
-
- // Initialize artifical cost
- Cost art_cost;
- if (std::numeric_limits<Cost>::is_exact) {
- art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
- } else {
- art_cost = std::numeric_limits<Cost>::min();
- for (int i = 0; i != _arc_num; ++i) {
- if (_cost[i] > art_cost) art_cost = _cost[i];
- }
- art_cost = (art_cost + 1) * _node_num;
+ ART_COST = (ART_COST + 1) * _node_num;
}
- // Run Circulation to check if a feasible solution exists
- typedef ConstMap<Arc, Flow> ConstArcMap;
- ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
- FlowNodeMap *csup = NULL;
- bool local_csup = false;
- if (_psupply) {
- csup = _psupply;
- } else {
- csup = new FlowNodeMap(_graph, 0);
- (*csup)[_psource] = _pstflow;
- (*csup)[_ptarget] = -_pstflow;
- local_csup = true;
+ // Initialize arc maps
+ for (int i = 0; i != _arc_num; ++i) {
+ _flow[i] = 0;
+ _state[i] = STATE_LOWER;
}
- bool circ_result = false;
- if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
- // GEQ problem type
- if (_plower) {
- if (_pupper) {
- Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
- circ(_graph, *_plower, *_pupper, *csup);
- circ_result = circ.run();
- } else {
- Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
- circ(_graph, *_plower, inf_arc_map, *csup);
- circ_result = circ.run();
- }
- } else {
- if (_pupper) {
- Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
- circ(_graph, zero_arc_map, *_pupper, *csup);
- circ_result = circ.run();
- } else {
- Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
- circ(_graph, zero_arc_map, inf_arc_map, *csup);
- circ_result = circ.run();
- }
- }
- } else {
- // LEQ problem type
- typedef ReverseDigraph<const GR> RevGraph;
- typedef NegMap<FlowNodeMap> NegNodeMap;
- RevGraph rgraph(_graph);
- NegNodeMap neg_csup(*csup);
- if (_plower) {
- if (_pupper) {
- Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
- circ(rgraph, *_plower, *_pupper, neg_csup);
- circ_result = circ.run();
- } else {
- Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
- circ(rgraph, *_plower, inf_arc_map, neg_csup);
- circ_result = circ.run();
- }
- } else {
- if (_pupper) {
- Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
- circ(rgraph, zero_arc_map, *_pupper, neg_csup);
- circ_result = circ.run();
- } else {
- Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
- circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
- circ_result = circ.run();
- }
- }
- }
- if (local_csup) delete csup;
- if (!circ_result) return false;
-
+
// Set data for the artificial root node
_root = _node_num;
_parent[_root] = -1;
_pred[_root] = -1;
_thread[_root] = 0;
_rev_thread[0] = _root;
- _succ_num[_root] = all_node_num;
+ _succ_num[_root] = _node_num + 1;
_last_succ[_root] = _root - 1;
- _supply[_root] = -sum_supply;
- if (sum_supply < 0) {
- _pi[_root] = -art_cost;
- } else {
- _pi[_root] = art_cost;
- }
-
- // Store the arcs in a mixed order
- int k = std::max(int(std::sqrt(double(_arc_num))), 10);
- int i = 0;
- for (ArcIt e(_graph); e != INVALID; ++e) {
- _arc_ref[i] = e;
- if ((i += k) >= _arc_num) i = (i % k) + 1;
- }
-
- // Initialize arc maps
- if (_pupper && _pcost) {
- for (int i = 0; i != _arc_num; ++i) {
- Arc e = _arc_ref[i];
- _source[i] = _node_id[_graph.source(e)];
- _target[i] = _node_id[_graph.target(e)];
- _cap[i] = (*_pupper)[e];
- _cost[i] = (*_pcost)[e];
- _flow[i] = 0;
- _state[i] = STATE_LOWER;
- }
- } else {
- for (int i = 0; i != _arc_num; ++i) {
- Arc e = _arc_ref[i];
- _source[i] = _node_id[_graph.source(e)];
- _target[i] = _node_id[_graph.target(e)];
- _flow[i] = 0;
- _state[i] = STATE_LOWER;
- }
- if (_pupper) {
- for (int i = 0; i != _arc_num; ++i)
- _cap[i] = (*_pupper)[_arc_ref[i]];
- } else {
- for (int i = 0; i != _arc_num; ++i)
- _cap[i] = inf_cap;
- }
- if (_pcost) {
- for (int i = 0; i != _arc_num; ++i)
- _cost[i] = (*_pcost)[_arc_ref[i]];
- } else {
- for (int i = 0; i != _arc_num; ++i)
- _cost[i] = 1;
- }
- }
-
- // Remove non-zero lower bounds
- if (_plower) {
- for (int i = 0; i != _arc_num; ++i) {
- Flow c = (*_plower)[_arc_ref[i]];
- if (c != 0) {
- _cap[i] -= c;
- _supply[_source[i]] -= c;
- _supply[_target[i]] += c;
- }
- }
- }
+ _supply[_root] = -_sum_supply;
+ _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
// Add artificial arcs and initialize the spanning tree data structure
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+ _parent[u] = _root;
+ _pred[u] = e;
_thread[u] = u + 1;
_rev_thread[u + 1] = u;
_succ_num[u] = 1;
_last_succ[u] = u;
- _parent[u] = _root;
- _pred[u] = e;
- _cost[e] = art_cost;
- _cap[e] = inf_cap;
+ _cost[e] = ART_COST;
+ _cap[e] = INF;
_state[e] = STATE_TREE;
- if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
+ if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
_flow[e] = _supply[u];
_forward[u] = true;
- _pi[u] = -art_cost + _pi[_root];
+ _pi[u] = -ART_COST + _pi[_root];
} else {
_flow[e] = -_supply[u];
_forward[u] = false;
- _pi[u] = art_cost + _pi[_root];
+ _pi[u] = ART_COST + _pi[_root];
}
}
@@ -1336,13 +1148,14 @@
}
delta = _cap[in_arc];
int result = 0;
- Flow d;
+ Value d;
int e;
// Search the cycle along the path form the first node to the root
for (int u = first; u != join; u = _parent[u]) {
e = _pred[u];
- d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
+ d = _forward[u] ?
+ _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
if (d < delta) {
delta = d;
u_out = u;
@@ -1352,7 +1165,8 @@
// Search the cycle along the path form the second node to the root
for (int u = second; u != join; u = _parent[u]) {
e = _pred[u];
- d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
+ d = _forward[u] ?
+ (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
if (d <= delta) {
delta = d;
u_out = u;
@@ -1374,7 +1188,7 @@
void changeFlow(bool change) {
// Augment along the cycle
if (delta > 0) {
- Flow val = _state[in_arc] * delta;
+ Value val = _state[in_arc] * delta;
_flow[in_arc] += val;
for (int u = _source[in_arc]; u != join; u = _parent[u]) {
_flow[_pred[u]] += _forward[u] ? -val : val;
@@ -1526,7 +1340,7 @@
}
// Execute the algorithm
- bool start(PivotRule pivot_rule) {
+ ProblemType start(PivotRule pivot_rule) {
// Select the pivot rule implementation
switch (pivot_rule) {
case FIRST_ELIGIBLE:
@@ -1540,41 +1354,55 @@
case ALTERING_LIST:
return start<AlteringListPivotRule>();
}
- return false;
+ return INFEASIBLE; // avoid warning
}
template <typename PivotRuleImpl>
- bool start() {
+ ProblemType start() {
PivotRuleImpl pivot(*this);
// Execute the Network Simplex algorithm
while (pivot.findEnteringArc()) {
findJoinNode();
bool change = findLeavingArc();
+ if (delta >= INF) return UNBOUNDED;
changeFlow(change);
if (change) {
updateTreeStructure();
updatePotential();
}
}
-
- // Copy flow values to _flow_map
- if (_plower) {
- for (int i = 0; i != _arc_num; ++i) {
- Arc e = _arc_ref[i];
- _flow_map->set(e, (*_plower)[e] + _flow[i]);
- }
- } else {
- for (int i = 0; i != _arc_num; ++i) {
- _flow_map->set(_arc_ref[i], _flow[i]);
+
+ // Check feasibility
+ if (_sum_supply < 0) {
+ for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+ if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
}
}
- // Copy potential values to _potential_map
- for (NodeIt n(_graph); n != INVALID; ++n) {
- _potential_map->set(n, _pi[_node_id[n]]);
+ else if (_sum_supply > 0) {
+ for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+ if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
+ }
+ }
+ else {
+ for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+ if (_flow[e] != 0) return INFEASIBLE;
+ }
}
- return true;
+ // Transform the solution and the supply map to the original form
+ if (_have_lower) {
+ for (int i = 0; i != _arc_num; ++i) {
+ Value c = _lower[i];
+ if (c != 0) {
+ _flow[i] += c;
+ _supply[_source[i]] += c;
+ _supply[_target[i]] -= c;
+ }
+ }
+ }
+
+ return OPTIMAL;
}
}; //class NetworkSimplex
diff -r cb8270a98660 -r e01957e96c67 lemon/preflow.h
--- a/lemon/preflow.h Wed Apr 29 16:15:29 2009 +0100
+++ b/lemon/preflow.h Wed Apr 29 19:22:14 2009 +0100
@@ -46,13 +46,13 @@
typedef CAP CapacityMap;
/// \brief The type of the flow values.
- typedef typename CapacityMap::Value Flow;
+ typedef typename CapacityMap::Value Value;
/// \brief The type of the map that stores the flow values.
///
/// The type of the map that stores the flow values.
/// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
- typedef typename Digraph::template ArcMap<Flow> FlowMap;
+ typedef typename Digraph::template ArcMap<Value> FlowMap;
/// \brief Instantiates a FlowMap.
///
@@ -84,7 +84,7 @@
/// \brief The tolerance used by the algorithm
///
/// The tolerance used by the algorithm to handle inexact computation.
- typedef lemon::Tolerance<Flow> Tolerance;
+ typedef lemon::Tolerance<Value> Tolerance;
};
@@ -125,7 +125,7 @@
///The type of the capacity map.
typedef typename Traits::CapacityMap CapacityMap;
///The type of the flow values.
- typedef typename Traits::Flow Flow;
+ typedef typename Traits::Value Value;
///The type of the flow map.
typedef typename Traits::FlowMap FlowMap;
@@ -151,7 +151,7 @@
Elevator* _level;
bool _local_level;
- typedef typename Digraph::template NodeMap<Flow> ExcessMap;
+ typedef typename Digraph::template NodeMap<Value> ExcessMap;
ExcessMap* _excess;
Tolerance _tolerance;
@@ -470,7 +470,7 @@
}
for (NodeIt n(_graph); n != INVALID; ++n) {
- Flow excess = 0;
+ Value excess = 0;
for (InArcIt e(_graph, n); e != INVALID; ++e) {
excess += (*_flow)[e];
}
@@ -519,7 +519,7 @@
_level->initFinish();
for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
- Flow rem = (*_capacity)[e] - (*_flow)[e];
+ Value rem = (*_capacity)[e] - (*_flow)[e];
if (_tolerance.positive(rem)) {
Node u = _graph.target(e);
if ((*_level)[u] == _level->maxLevel()) continue;
@@ -531,7 +531,7 @@
}
}
for (InArcIt e(_graph, _source); e != INVALID; ++e) {
- Flow rem = (*_flow)[e];
+ Value rem = (*_flow)[e];
if (_tolerance.positive(rem)) {
Node v = _graph.source(e);
if ((*_level)[v] == _level->maxLevel()) continue;
@@ -564,11 +564,11 @@
int num = _node_num;
while (num > 0 && n != INVALID) {
- Flow excess = (*_excess)[n];
+ Value excess = (*_excess)[n];
int new_level = _level->maxLevel();
for (OutArcIt e(_graph, n); e != INVALID; ++e) {
- Flow rem = (*_capacity)[e] - (*_flow)[e];
+ Value rem = (*_capacity)[e] - (*_flow)[e];
if (!_tolerance.positive(rem)) continue;
Node v = _graph.target(e);
if ((*_level)[v] < level) {
@@ -591,7 +591,7 @@
}
for (InArcIt e(_graph, n); e != INVALID; ++e) {
- Flow rem = (*_flow)[e];
+ Value rem = (*_flow)[e];
if (!_tolerance.positive(rem)) continue;
Node v = _graph.source(e);
if ((*_level)[v] < level) {
@@ -637,11 +637,11 @@
num = _node_num * 20;
while (num > 0 && n != INVALID) {
- Flow excess = (*_excess)[n];
+ Value excess = (*_excess)[n];
int new_level = _level->maxLevel();
for (OutArcIt e(_graph, n); e != INVALID; ++e) {
- Flow rem = (*_capacity)[e] - (*_flow)[e];
+ Value rem = (*_capacity)[e] - (*_flow)[e];
if (!_tolerance.positive(rem)) continue;
Node v = _graph.target(e);
if ((*_level)[v] < level) {
@@ -664,7 +664,7 @@
}
for (InArcIt e(_graph, n); e != INVALID; ++e) {
- Flow rem = (*_flow)[e];
+ Value rem = (*_flow)[e];
if (!_tolerance.positive(rem)) continue;
Node v = _graph.source(e);
if ((*_level)[v] < level) {
@@ -778,12 +778,12 @@
Node n;
while ((n = _level->highestActive()) != INVALID) {
- Flow excess = (*_excess)[n];
+ Value excess = (*_excess)[n];
int level = _level->highestActiveLevel();
int new_level = _level->maxLevel();
for (OutArcIt e(_graph, n); e != INVALID; ++e) {
- Flow rem = (*_capacity)[e] - (*_flow)[e];
+ Value rem = (*_capacity)[e] - (*_flow)[e];
if (!_tolerance.positive(rem)) continue;
Node v = _graph.target(e);
if ((*_level)[v] < level) {
@@ -806,7 +806,7 @@
}
for (InArcIt e(_graph, n); e != INVALID; ++e) {
- Flow rem = (*_flow)[e];
+ Value rem = (*_flow)[e];
if (!_tolerance.positive(rem)) continue;
Node v = _graph.source(e);
if ((*_level)[v] < level) {
@@ -897,18 +897,18 @@
///
/// \pre Either \ref run() or \ref init() must be called before
/// using this function.
- Flow flowValue() const {
+ Value flowValue() const {
return (*_excess)[_target];
}
- /// \brief Returns the flow on the given arc.
+ /// \brief Returns the flow value on the given arc.
///
- /// Returns the flow on the given arc. This method can
+ /// Returns the flow value on the given arc. This method can
/// be called after the second phase of the algorithm.
///
/// \pre Either \ref run() or \ref init() must be called before
/// using this function.
- Flow flow(const Arc& arc) const {
+ Value flow(const Arc& arc) const {
return (*_flow)[arc];
}
diff -r cb8270a98660 -r e01957e96c67 test/lp_test.cc
--- a/test/lp_test.cc Wed Apr 29 16:15:29 2009 +0100
+++ b/test/lp_test.cc Wed Apr 29 19:22:14 2009 +0100
@@ -21,9 +21,7 @@
#include "test_tools.h"
#include <lemon/tolerance.h>
-#ifdef HAVE_CONFIG_H
#include <lemon/config.h>
-#endif
#ifdef LEMON_HAVE_GLPK
#include <lemon/glpk.h>
diff -r cb8270a98660 -r e01957e96c67 test/min_cost_flow_test.cc
--- a/test/min_cost_flow_test.cc Wed Apr 29 16:15:29 2009 +0100
+++ b/test/min_cost_flow_test.cc Wed Apr 29 19:22:14 2009 +0100
@@ -18,6 +18,7 @@
#include <iostream>
#include <fstream>
+#include <limits>
#include <lemon/list_graph.h>
#include <lemon/lgf_reader.h>
@@ -33,57 +34,57 @@
char test_lgf[] =
"@nodes\n"
- "label sup1 sup2 sup3 sup4 sup5\n"
- " 1 20 27 0 20 30\n"
- " 2 -4 0 0 -8 -3\n"
- " 3 0 0 0 0 0\n"
- " 4 0 0 0 0 0\n"
- " 5 9 0 0 6 11\n"
- " 6 -6 0 0 -5 -6\n"
- " 7 0 0 0 0 0\n"
- " 8 0 0 0 0 3\n"
- " 9 3 0 0 0 0\n"
- " 10 -2 0 0 -7 -2\n"
- " 11 0 0 0 -10 0\n"
- " 12 -20 -27 0 -30 -20\n"
- "\n"
+ "label sup1 sup2 sup3 sup4 sup5 sup6\n"
+ " 1 20 27 0 30 20 30\n"
+ " 2 -4 0 0 0 -8 -3\n"
+ " 3 0 0 0 0 0 0\n"
+ " 4 0 0 0 0 0 0\n"
+ " 5 9 0 0 0 6 11\n"
+ " 6 -6 0 0 0 -5 -6\n"
+ " 7 0 0 0 0 0 0\n"
+ " 8 0 0 0 0 0 3\n"
+ " 9 3 0 0 0 0 0\n"
+ " 10 -2 0 0 0 -7 -2\n"
+ " 11 0 0 0 0 -10 0\n"
+ " 12 -20 -27 0 -30 -30 -20\n"
+ "\n"
"@arcs\n"
- " cost cap low1 low2\n"
- " 1 2 70 11 0 8\n"
- " 1 3 150 3 0 1\n"
- " 1 4 80 15 0 2\n"
- " 2 8 80 12 0 0\n"
- " 3 5 140 5 0 3\n"
- " 4 6 60 10 0 1\n"
- " 4 7 80 2 0 0\n"
- " 4 8 110 3 0 0\n"
- " 5 7 60 14 0 0\n"
- " 5 11 120 12 0 0\n"
- " 6 3 0 3 0 0\n"
- " 6 9 140 4 0 0\n"
- " 6 10 90 8 0 0\n"
- " 7 1 30 5 0 0\n"
- " 8 12 60 16 0 4\n"
- " 9 12 50 6 0 0\n"
- "10 12 70 13 0 5\n"
- "10 2 100 7 0 0\n"
- "10 7 60 10 0 0\n"
- "11 10 20 14 0 6\n"
- "12 11 30 10 0 0\n"
+ " cost cap low1 low2 low3\n"
+ " 1 2 70 11 0 8 8\n"
+ " 1 3 150 3 0 1 0\n"
+ " 1 4 80 15 0 2 2\n"
+ " 2 8 80 12 0 0 0\n"
+ " 3 5 140 5 0 3 1\n"
+ " 4 6 60 10 0 1 0\n"
+ " 4 7 80 2 0 0 0\n"
+ " 4 8 110 3 0 0 0\n"
+ " 5 7 60 14 0 0 0\n"
+ " 5 11 120 12 0 0 0\n"
+ " 6 3 0 3 0 0 0\n"
+ " 6 9 140 4 0 0 0\n"
+ " 6 10 90 8 0 0 0\n"
+ " 7 1 30 5 0 0 -5\n"
+ " 8 12 60 16 0 4 3\n"
+ " 9 12 50 6 0 0 0\n"
+ "10 12 70 13 0 5 2\n"
+ "10 2 100 7 0 0 0\n"
+ "10 7 60 10 0 0 -3\n"
+ "11 10 20 14 0 6 -20\n"
+ "12 11 30 10 0 0 -10\n"
"\n"
"@attributes\n"
"source 1\n"
"target 12\n";
-enum ProblemType {
+enum SupplyType {
EQ,
GEQ,
LEQ
};
// Check the interface of an MCF algorithm
-template <typename GR, typename Flow, typename Cost>
+template <typename GR, typename Value, typename Cost>
class McfClassConcept
{
public:
@@ -94,53 +95,46 @@
checkConcept<concepts::Digraph, GR>();
MCF mcf(g);
+ const MCF& const_mcf = mcf;
b = mcf.reset()
.lowerMap(lower)
.upperMap(upper)
- .capacityMap(upper)
- .boundMaps(lower, upper)
.costMap(cost)
.supplyMap(sup)
.stSupply(n, n, k)
- .flowMap(flow)
- .potentialMap(pot)
.run();
-
- const MCF& const_mcf = mcf;
- const typename MCF::FlowMap &fm = const_mcf.flowMap();
- const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
-
- v = const_mcf.totalCost();
- double x = const_mcf.template totalCost<double>();
+ c = const_mcf.totalCost();
+ x = const_mcf.template totalCost<double>();
v = const_mcf.flow(a);
- v = const_mcf.potential(n);
-
- ignore_unused_variable_warning(fm);
- ignore_unused_variable_warning(pm);
- ignore_unused_variable_warning(x);
+ c = const_mcf.potential(n);
+ const_mcf.flowMap(fm);
+ const_mcf.potentialMap(pm);
}
typedef typename GR::Node Node;
typedef typename GR::Arc Arc;
- typedef concepts::ReadMap<Node, Flow> NM;
- typedef concepts::ReadMap<Arc, Flow> FAM;
+ typedef concepts::ReadMap<Node, Value> NM;
+ typedef concepts::ReadMap<Arc, Value> VAM;
typedef concepts::ReadMap<Arc, Cost> CAM;
+ typedef concepts::WriteMap<Arc, Value> FlowMap;
+ typedef concepts::WriteMap<Node, Cost> PotMap;
const GR &g;
- const FAM &lower;
- const FAM &upper;
+ const VAM &lower;
+ const VAM &upper;
const CAM &cost;
const NM ⊃
const Node &n;
const Arc &a;
- const Flow &k;
- Flow v;
+ const Value &k;
+ FlowMap fm;
+ PotMap pm;
bool b;
-
- typename MCF::FlowMap &flow;
- typename MCF::PotentialMap &pot;
+ double x;
+ typename MCF::Value v;
+ typename MCF::Cost c;
};
};
@@ -151,7 +145,7 @@
typename SM, typename FM >
bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
const SM& supply, const FM& flow,
- ProblemType type = EQ )
+ SupplyType type = EQ )
{
TEMPLATE_DIGRAPH_TYPEDEFS(GR);
@@ -208,21 +202,25 @@
// Run a minimum cost flow algorithm and check the results
template < typename MCF, typename GR,
typename LM, typename UM,
- typename CM, typename SM >
-void checkMcf( const MCF& mcf, bool mcf_result,
+ typename CM, typename SM,
+ typename PT >
+void checkMcf( const MCF& mcf, PT mcf_result,
const GR& gr, const LM& lower, const UM& upper,
const CM& cost, const SM& supply,
- bool result, typename CM::Value total,
+ PT result, bool optimal, typename CM::Value total,
const std::string &test_id = "",
- ProblemType type = EQ )
+ SupplyType type = EQ )
{
check(mcf_result == result, "Wrong result " + test_id);
- if (result) {
- check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
+ if (optimal) {
+ typename GR::template ArcMap<typename SM::Value> flow(gr);
+ typename GR::template NodeMap<typename CM::Value> pi(gr);
+ mcf.flowMap(flow);
+ mcf.potentialMap(pi);
+ check(checkFlow(gr, lower, upper, supply, flow, type),
"The flow is not feasible " + test_id);
check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
- check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
- mcf.potentialMap()),
+ check(checkPotential(gr, lower, upper, cost, supply, flow, pi),
"Wrong potentials " + test_id);
}
}
@@ -231,11 +229,13 @@
{
// Check the interfaces
{
- typedef int Flow;
- typedef int Cost;
typedef concepts::Digraph GR;
- checkConcept< McfClassConcept<GR, Flow, Cost>,
- NetworkSimplex<GR, Flow, Cost> >();
+ checkConcept< McfClassConcept<GR, int, int>,
+ NetworkSimplex<GR> >();
+ checkConcept< McfClassConcept<GR, double, double>,
+ NetworkSimplex<GR, double> >();
+ checkConcept< McfClassConcept<GR, int, double>,
+ NetworkSimplex<GR, int, double> >();
}
// Run various MCF tests
@@ -244,8 +244,8 @@
// Read the test digraph
Digraph gr;
- Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
- Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
+ Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
+ Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
Node v, w;
@@ -255,14 +255,56 @@
.arcMap("cap", u)
.arcMap("low1", l1)
.arcMap("low2", l2)
+ .arcMap("low3", l3)
.nodeMap("sup1", s1)
.nodeMap("sup2", s2)
.nodeMap("sup3", s3)
.nodeMap("sup4", s4)
.nodeMap("sup5", s5)
+ .nodeMap("sup6", s6)
.node("source", v)
.node("target", w)
.run();
+
+ // Build a test digraph for testing negative costs
+ Digraph ngr;
+ Node n1 = ngr.addNode();
+ Node n2 = ngr.addNode();
+ Node n3 = ngr.addNode();
+ Node n4 = ngr.addNode();
+ Node n5 = ngr.addNode();
+ Node n6 = ngr.addNode();
+ Node n7 = ngr.addNode();
+
+ Arc a1 = ngr.addArc(n1, n2);
+ Arc a2 = ngr.addArc(n1, n3);
+ Arc a3 = ngr.addArc(n2, n4);
+ Arc a4 = ngr.addArc(n3, n4);
+ Arc a5 = ngr.addArc(n3, n2);
+ Arc a6 = ngr.addArc(n5, n3);
+ Arc a7 = ngr.addArc(n5, n6);
+ Arc a8 = ngr.addArc(n6, n7);
+ Arc a9 = ngr.addArc(n7, n5);
+
+ Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
+ ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
+ Digraph::NodeMap<int> ns(ngr, 0);
+
+ nl2[a7] = 1000;
+ nl2[a8] = -1000;
+
+ ns[n1] = 100;
+ ns[n4] = -100;
+
+ nc[a1] = 100;
+ nc[a2] = 30;
+ nc[a3] = 20;
+ nc[a4] = 80;
+ nc[a5] = 50;
+ nc[a6] = 10;
+ nc[a7] = 80;
+ nc[a8] = 30;
+ nc[a9] = -120;
// A. Test NetworkSimplex with the default pivot rule
{
@@ -271,63 +313,77 @@
// Check the equality form
mcf.upperMap(u).costMap(c);
checkMcf(mcf, mcf.supplyMap(s1).run(),
- gr, l1, u, c, s1, true, 5240, "#A1");
+ gr, l1, u, c, s1, mcf.OPTIMAL, true, 5240, "#A1");
checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
- gr, l1, u, c, s2, true, 7620, "#A2");
+ gr, l1, u, c, s2, mcf.OPTIMAL, true, 7620, "#A2");
mcf.lowerMap(l2);
checkMcf(mcf, mcf.supplyMap(s1).run(),
- gr, l2, u, c, s1, true, 5970, "#A3");
+ gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#A3");
checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
- gr, l2, u, c, s2, true, 8010, "#A4");
+ gr, l2, u, c, s2, mcf.OPTIMAL, true, 8010, "#A4");
mcf.reset();
checkMcf(mcf, mcf.supplyMap(s1).run(),
- gr, l1, cu, cc, s1, true, 74, "#A5");
+ gr, l1, cu, cc, s1, mcf.OPTIMAL, true, 74, "#A5");
checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
- gr, l2, cu, cc, s2, true, 94, "#A6");
+ gr, l2, cu, cc, s2, mcf.OPTIMAL, true, 94, "#A6");
mcf.reset();
checkMcf(mcf, mcf.run(),
- gr, l1, cu, cc, s3, true, 0, "#A7");
- checkMcf(mcf, mcf.boundMaps(l2, u).run(),
- gr, l2, u, cc, s3, false, 0, "#A8");
+ gr, l1, cu, cc, s3, mcf.OPTIMAL, true, 0, "#A7");
+ checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
+ gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
+ mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
+ checkMcf(mcf, mcf.run(),
+ gr, l3, u, c, s4, mcf.OPTIMAL, true, 6360, "#A9");
// Check the GEQ form
- mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
+ mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
checkMcf(mcf, mcf.run(),
- gr, l1, u, c, s4, true, 3530, "#A9", GEQ);
- mcf.problemType(mcf.GEQ);
+ gr, l1, u, c, s5, mcf.OPTIMAL, true, 3530, "#A10", GEQ);
+ mcf.supplyType(mcf.GEQ);
checkMcf(mcf, mcf.lowerMap(l2).run(),
- gr, l2, u, c, s4, true, 4540, "#A10", GEQ);
- mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
+ gr, l2, u, c, s5, mcf.OPTIMAL, true, 4540, "#A11", GEQ);
+ mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
checkMcf(mcf, mcf.run(),
- gr, l2, u, c, s5, false, 0, "#A11", GEQ);
+ gr, l2, u, c, s6, mcf.INFEASIBLE, false, 0, "#A12", GEQ);
// Check the LEQ form
- mcf.reset().problemType(mcf.LEQ);
- mcf.upperMap(u).costMap(c).supplyMap(s5);
+ mcf.reset().supplyType(mcf.LEQ);
+ mcf.upperMap(u).costMap(c).supplyMap(s6);
checkMcf(mcf, mcf.run(),
- gr, l1, u, c, s5, true, 5080, "#A12", LEQ);
+ gr, l1, u, c, s6, mcf.OPTIMAL, true, 5080, "#A13", LEQ);
checkMcf(mcf, mcf.lowerMap(l2).run(),
- gr, l2, u, c, s5, true, 5930, "#A13", LEQ);
- mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
+ gr, l2, u, c, s6, mcf.OPTIMAL, true, 5930, "#A14", LEQ);
+ mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
checkMcf(mcf, mcf.run(),
- gr, l2, u, c, s4, false, 0, "#A14", LEQ);
+ gr, l2, u, c, s5, mcf.INFEASIBLE, false, 0, "#A15", LEQ);
+
+ // Check negative costs
+ NetworkSimplex<Digraph> nmcf(ngr);
+ nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
+ checkMcf(nmcf, nmcf.run(),
+ ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false, 0, "#A16");
+ checkMcf(nmcf, nmcf.upperMap(nu2).run(),
+ ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true, -40000, "#A17");
+ nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
+ checkMcf(nmcf, nmcf.run(),
+ ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false, 0, "#A18");
}
// B. Test NetworkSimplex with each pivot rule
{
NetworkSimplex<Digraph> mcf(gr);
- mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
+ mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
- gr, l2, u, c, s1, true, 5970, "#B1");
+ gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B1");
checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
- gr, l2, u, c, s1, true, 5970, "#B2");
+ gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B2");
checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
- gr, l2, u, c, s1, true, 5970, "#B3");
+ gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B3");
checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
- gr, l2, u, c, s1, true, 5970, "#B4");
+ gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B4");
checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
- gr, l2, u, c, s1, true, 5970, "#B5");
+ gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B5");
}
return 0;
diff -r cb8270a98660 -r e01957e96c67 test/mip_test.cc
--- a/test/mip_test.cc Wed Apr 29 16:15:29 2009 +0100
+++ b/test/mip_test.cc Wed Apr 29 19:22:14 2009 +0100
@@ -18,9 +18,7 @@
#include "test_tools.h"
-#ifdef HAVE_CONFIG_H
#include <lemon/config.h>
-#endif
#ifdef LEMON_HAVE_CPLEX
#include <lemon/cplex.h>
diff -r cb8270a98660 -r e01957e96c67 tools/dimacs-solver.cc
--- a/tools/dimacs-solver.cc Wed Apr 29 16:15:29 2009 +0100
+++ b/tools/dimacs-solver.cc Wed Apr 29 19:22:14 2009 +0100
@@ -119,8 +119,8 @@
ti.restart();
NetworkSimplex<Digraph, Value> ns(g);
- ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
- if (sum_sup > 0) ns.problemType(ns.LEQ);
+ ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup);
+ if (sum_sup > 0) ns.supplyType(ns.LEQ);
if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';
ti.restart();
bool res = ns.run();