↑ Collapse diff ↑
Show white space 6 line context
... ...
@@ -17,8 +17,6 @@
17 17
FIND_PACKAGE(CPLEX)
18 18
FIND_PACKAGE(COIN)
19 19

	
20
ADD_DEFINITIONS(-DHAVE_CONFIG_H)
21

	
22 20
IF(MSVC)
23 21
  SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4250 /wd4355 /wd4800 /wd4996")
24 22
# Suppressed warnings:
... ...
@@ -28,8 +26,6 @@
28 26
# C4996: 'function': was declared deprecated
29 27
ENDIF(MSVC)
30 28

	
31
ADD_DEFINITIONS(-DHAVE_CONFIG_H)
32

	
33 29
INCLUDE(CheckTypeSize)
34 30
CHECK_TYPE_SIZE("long long" LEMON_LONG_LONG)
35 31

	
Show white space 6 line context
... ...
@@ -11,11 +11,12 @@
11 11
	m4/lx_check_cplex.m4 \
12 12
	m4/lx_check_glpk.m4 \
13 13
	m4/lx_check_soplex.m4 \
14
	m4/lx_check_clp.m4 \
15
	m4/lx_check_cbc.m4 \
14
	m4/lx_check_coin.m4 \
16 15
	CMakeLists.txt \
17 16
	cmake/FindGhostscript.cmake \
17
	cmake/FindCPLEX.cmake \
18 18
	cmake/FindGLPK.cmake \
19
	cmake/FindCOIN.cmake \
19 20
	cmake/version.cmake.in \
20 21
	cmake/version.cmake \
21 22
	cmake/nsis/lemon.ico \
Show white space 6 line context
1 1
SET(COIN_ROOT_DIR "" CACHE PATH "COIN root directory")
2 2

	
3 3
FIND_PATH(COIN_INCLUDE_DIR coin/CoinUtilsConfig.h
4
  PATHS ${COIN_ROOT_DIR}/include)
5

	
6
FIND_LIBRARY(COIN_CBC_LIBRARY libCbc
7
  PATHS ${COIN_ROOT_DIR}/lib)
8
FIND_LIBRARY(COIN_CBC_SOLVER_LIBRARY libCbcSolver
9
  PATHS ${COIN_ROOT_DIR}/lib)
10
FIND_LIBRARY(COIN_CGL_LIBRARY libCgl
11
  PATHS ${COIN_ROOT_DIR}/lib)
12
FIND_LIBRARY(COIN_CLP_LIBRARY libClp
13
  PATHS ${COIN_ROOT_DIR}/lib)
14
FIND_LIBRARY(COIN_COIN_UTILS_LIBRARY libCoinUtils
15
  PATHS ${COIN_ROOT_DIR}/lib)
16
FIND_LIBRARY(COIN_OSI_LIBRARY libOsi
17
  PATHS ${COIN_ROOT_DIR}/lib)
18
FIND_LIBRARY(COIN_OSI_CBC_LIBRARY libOsiCbc
19
  PATHS ${COIN_ROOT_DIR}/lib)
20
FIND_LIBRARY(COIN_OSI_CLP_LIBRARY libOsiClp
21
  PATHS ${COIN_ROOT_DIR}/lib)
22
FIND_LIBRARY(COIN_OSI_VOL_LIBRARY libOsiVol
23
  PATHS ${COIN_ROOT_DIR}/lib)
24
FIND_LIBRARY(COIN_VOL_LIBRARY libVol
25
  PATHS ${COIN_ROOT_DIR}/lib)
4
  HINTS ${COIN_ROOT_DIR}/include
5
)
6
FIND_LIBRARY(COIN_CBC_LIBRARY
7
  NAMES Cbc libCbc
8
  HINTS ${COIN_ROOT_DIR}/lib
9
)
10
FIND_LIBRARY(COIN_CBC_SOLVER_LIBRARY
11
  NAMES CbcSolver libCbcSolver
12
  HINTS ${COIN_ROOT_DIR}/lib
13
)
14
FIND_LIBRARY(COIN_CGL_LIBRARY
15
  NAMES Cgl libCgl
16
  HINTS ${COIN_ROOT_DIR}/lib
17
)
18
FIND_LIBRARY(COIN_CLP_LIBRARY
19
  NAMES Clp libClp
20
  HINTS ${COIN_ROOT_DIR}/lib
21
)
22
FIND_LIBRARY(COIN_COIN_UTILS_LIBRARY
23
  NAMES CoinUtils libCoinUtils
24
  HINTS ${COIN_ROOT_DIR}/lib
25
)
26
FIND_LIBRARY(COIN_OSI_LIBRARY
27
  NAMES Osi libOsi
28
  HINTS ${COIN_ROOT_DIR}/lib
29
)
30
FIND_LIBRARY(COIN_OSI_CBC_LIBRARY
31
  NAMES OsiCbc libOsiCbc
32
  HINTS ${COIN_ROOT_DIR}/lib
33
)
34
FIND_LIBRARY(COIN_OSI_CLP_LIBRARY
35
  NAMES OsiClp libOsiClp
36
  HINTS ${COIN_ROOT_DIR}/lib
37
)
38
FIND_LIBRARY(COIN_OSI_VOL_LIBRARY
39
  NAMES OsiVol libOsiVol
40
  HINTS ${COIN_ROOT_DIR}/lib
41
)
42
FIND_LIBRARY(COIN_VOL_LIBRARY
43
  NAMES Vol libVol
44
  HINTS ${COIN_ROOT_DIR}/lib
45
)
26 46

	
27 47
INCLUDE(FindPackageHandleStandardArgs)
28 48
FIND_PACKAGE_HANDLE_STANDARD_ARGS(COIN DEFAULT_MSG
Show white space 6 line context
1
SET(CPLEX_ROOT_DIR "" CACHE PATH "CPLEX root directory")
2

	
1 3
FIND_PATH(CPLEX_INCLUDE_DIR
2 4
  ilcplex/cplex.h
3
  PATHS "C:/ILOG/CPLEX91/include")
4

	
5
  PATHS "C:/ILOG/CPLEX91/include"
6
  PATHS "/opt/ilog/cplex91/include"
7
  HINTS ${CPLEX_ROOT_DIR}/include
8
)
5 9
FIND_LIBRARY(CPLEX_LIBRARY
6
  NAMES cplex91
7
  PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda")
10
  cplex91
11
  PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda"
12
  PATHS "/opt/ilog/cplex91/bin"
13
  HINTS ${CPLEX_ROOT_DIR}/bin
14
)
8 15

	
9 16
INCLUDE(FindPackageHandleStandardArgs)
10 17
FIND_PACKAGE_HANDLE_STANDARD_ARGS(CPLEX DEFAULT_MSG CPLEX_LIBRARY CPLEX_INCLUDE_DIR)
11 18

	
12 19
FIND_PATH(CPLEX_BIN_DIR
13 20
  cplex91.dll
14
  PATHS "C:/ILOG/CPLEX91/bin/x86_win32")
21
  PATHS "C:/ILOG/CPLEX91/bin/x86_win32"
22
)
15 23

	
16 24
IF(CPLEX_FOUND)
17 25
  SET(CPLEX_INCLUDE_DIRS ${CPLEX_INCLUDE_DIR})
18 26
  SET(CPLEX_LIBRARIES ${CPLEX_LIBRARY})
27
  IF(CMAKE_SYSTEM_NAME STREQUAL "Linux")
28
    SET(CPLEX_LIBRARIES "${CPLEX_LIBRARIES};m;pthread")
29
  ENDIF(CMAKE_SYSTEM_NAME STREQUAL "Linux")
19 30
ENDIF(CPLEX_FOUND)
20 31

	
21 32
MARK_AS_ADVANCED(CPLEX_LIBRARY CPLEX_INCLUDE_DIR CPLEX_BIN_DIR)
Show white space 6 line context
1
SET(GLPK_ROOT_DIR "" CACHE PATH "GLPK root directory")
2

	
1 3
SET(GLPK_REGKEY "[HKEY_LOCAL_MACHINE\\SOFTWARE\\GnuWin32\\Glpk;InstallPath]")
2 4
GET_FILENAME_COMPONENT(GLPK_ROOT_PATH ${GLPK_REGKEY} ABSOLUTE)
3 5

	
4 6
FIND_PATH(GLPK_INCLUDE_DIR
5 7
  glpk.h
6
  PATHS ${GLPK_REGKEY}/include)
8
  PATHS ${GLPK_REGKEY}/include
9
  HINTS ${GLPK_ROOT_DIR}/include
10
)
11
FIND_LIBRARY(GLPK_LIBRARY
12
  glpk
13
  PATHS ${GLPK_REGKEY}/lib
14
  HINTS ${GLPK_ROOT_DIR}/lib
15
)
7 16

	
8
FIND_LIBRARY(GLPK_LIBRARY
9
  NAMES glpk
10
  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)
16 50
  SET(GLPK_INCLUDE_DIRS ${GLPK_INCLUDE_DIR})
Show white space 6 line context
... ...
@@ -352,17 +352,17 @@
352 352
minimum total cost from a set of supply nodes to a set of demand nodes
353 353
in a network with capacity constraints (lower and upper bounds)
354 354
and arc costs.
355
Formally, let \f$G=(V,A)\f$ be a digraph,
356
\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
355
Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
356
\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
357 357
upper bounds for the flow values on the arcs, for which
358
\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
359
\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
360
on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
358
\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
359
\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
360
on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
361 361
signed supply values of the nodes.
362 362
If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
363 363
supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
364 364
\f$-sup(u)\f$ demand.
365
A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
365
A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
366 366
of the following optimization problem.
367 367

	
368 368
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
... ...
@@ -404,7 +404,7 @@
404 404

	
405 405
The dual solution of the minimum cost flow problem is represented by node 
406 406
potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
407
An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
407
An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
408 408
is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
409 409
node potentials the following \e complementary \e slackness optimality
410 410
conditions hold.
... ...
@@ -413,15 +413,15 @@
413 413
   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
414 414
   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
415 415
   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
416
 - For all \f$u\in V\f$:
416
 - For all \f$u\in V\f$ nodes:
417 417
   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
418 418
     then \f$\pi(u)=0\f$.
419 419
 
420 420
Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
421
\f$uv\in A\f$ with respect to the 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

	
427 427
LEMON contains several algorithms for solving minimum cost flow problems.
Show white space 6 line context
... ...
@@ -15,6 +15,7 @@
15 15
	lemon/random.cc \
16 16
	lemon/bits/windows.cc
17 17

	
18
nodist_lemon_HEADERS = lemon/config.h	
18 19

	
19 20
lemon_libemon_la_CXXFLAGS = \
20 21
	$(AM_CXXFLAGS) \
... ...
@@ -57,11 +58,11 @@
57 58
	lemon/assert.h \
58 59
	lemon/bfs.h \
59 60
	lemon/bin_heap.h \
61
	lemon/cbc.h \
60 62
	lemon/circulation.h \
61 63
	lemon/clp.h \
62 64
	lemon/color.h \
63 65
	lemon/concept_check.h \
64
	lemon/config.h \
65 66
	lemon/connectivity.h \
66 67
	lemon/counter.h \
67 68
	lemon/core.h \
Show white space 6 line context
... ...
@@ -64,15 +64,15 @@
64 64
    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
65 65
    typedef SM SupplyMap;
66 66

	
67
    /// \brief The type of the flow values.
68
    typedef typename SupplyMap::Value Flow;
67
    /// \brief The type of the flow and supply values.
68
    typedef typename SupplyMap::Value Value;
69 69

	
70 70
    /// \brief The type of the map that stores the flow values.
71 71
    ///
72 72
    /// The type of the map that stores the flow values.
73 73
    /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
74 74
    /// concept.
75
    typedef typename Digraph::template ArcMap<Flow> FlowMap;
75
    typedef typename Digraph::template ArcMap<Value> FlowMap;
76 76

	
77 77
    /// \brief Instantiates a FlowMap.
78 78
    ///
... ...
@@ -104,7 +104,7 @@
104 104
    /// \brief The tolerance used by the algorithm
105 105
    ///
106 106
    /// The tolerance used by the algorithm to handle inexact computation.
107
    typedef lemon::Tolerance<Flow> Tolerance;
107
    typedef lemon::Tolerance<Value> Tolerance;
108 108

	
109 109
  };
110 110

	
... ...
@@ -187,8 +187,8 @@
187 187
    typedef TR Traits;
188 188
    ///The type of the digraph the algorithm runs on.
189 189
    typedef typename Traits::Digraph Digraph;
190
    ///The type of the flow values.
191
    typedef typename Traits::Flow Flow;
190
    ///The type of the flow and supply values.
191
    typedef typename Traits::Value Value;
192 192

	
193 193
    ///The type of the lower bound map.
194 194
    typedef typename Traits::LowerMap LowerMap;
... ...
@@ -221,7 +221,7 @@
221 221
    Elevator* _level;
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

	
227 227
    Tolerance _tol;
... ...
@@ -530,7 +530,7 @@
530 530
          (*_excess)[_g.target(e)] += (*_lo)[e];
531 531
          (*_excess)[_g.source(e)] -= (*_lo)[e];
532 532
        } else {
533
          Flow fc = -(*_excess)[_g.target(e)];
533
          Value fc = -(*_excess)[_g.target(e)];
534 534
          _flow->set(e, fc);
535 535
          (*_excess)[_g.target(e)] = 0;
536 536
          (*_excess)[_g.source(e)] -= fc;
... ...
@@ -563,11 +563,11 @@
563 563
      while((act=_level->highestActive())!=INVALID) {
564 564
        int actlevel=(*_level)[act];
565 565
        int mlevel=_node_num;
566
        Flow exc=(*_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
          Flow fc=(*_up)[e]-(*_flow)[e];
570
          Value fc=(*_up)[e]-(*_flow)[e];
571 571
          if(!_tol.positive(fc)) continue;
572 572
          if((*_level)[v]<actlevel) {
573 573
            if(!_tol.less(fc, exc)) {
... ...
@@ -591,7 +591,7 @@
591 591
        }
592 592
        for(InArcIt e(_g,act);e!=INVALID; ++e) {
593 593
          Node v = _g.source(e);
594
          Flow fc=(*_flow)[e]-(*_lo)[e];
594
          Value fc=(*_flow)[e]-(*_lo)[e];
595 595
          if(!_tol.positive(fc)) continue;
596 596
          if((*_level)[v]<actlevel) {
597 597
            if(!_tol.less(fc, exc)) {
... ...
@@ -661,13 +661,13 @@
661 661

	
662 662
    ///@{
663 663

	
664
    /// \brief Returns the flow on the given arc.
664
    /// \brief Returns the flow value on the given arc.
665 665
    ///
666
    /// Returns the flow on the given arc.
666
    /// Returns the flow value on the given arc.
667 667
    ///
668 668
    /// \pre Either \ref run() or \ref init() must be called before
669 669
    /// using this function.
670
    Flow flow(const Arc& arc) const {
670
    Value flow(const Arc& arc) const {
671 671
      return (*_flow)[arc];
672 672
    }
673 673

	
... ...
@@ -750,7 +750,7 @@
750 750
        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
751 751
      for(NodeIt n(_g);n!=INVALID;++n)
752 752
        {
753
          Flow dif=-(*_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];
756 756
          if(_tol.negative(dif)) return false;
... ...
@@ -765,10 +765,10 @@
765 765
    ///\sa barrierMap()
766 766
    bool checkBarrier() const
767 767
    {
768
      Flow delta=0;
769
      Flow inf_cap = std::numeric_limits<Flow>::has_infinity ?
770
        std::numeric_limits<Flow>::infinity() :
771
        std::numeric_limits<Flow>::max();
768
      Value delta=0;
769
      Value inf_cap = std::numeric_limits<Value>::has_infinity ?
770
        std::numeric_limits<Value>::infinity() :
771
        std::numeric_limits<Value>::max();
772 772
      for(NodeIt n(_g);n!=INVALID;++n)
773 773
        if(barrier(n))
774 774
          delta-=(*_supply)[n];
Show white space 6 line context
... ...
@@ -22,7 +22,7 @@
22 22
#include <vector>
23 23
#include <algorithm>
24 24

	
25
#include <lemon/core.h>
25
#include <lemon/config.h>
26 26
#include <lemon/bits/enable_if.h>
27 27
#include <lemon/bits/traits.h>
28 28
#include <lemon/assert.h>
Show white space 6 line context
... ...
@@ -30,9 +30,6 @@
30 30

	
31 31
#include <lemon/core.h>
32 32
#include <lemon/math.h>
33
#include <lemon/maps.h>
34
#include <lemon/circulation.h>
35
#include <lemon/adaptors.h>
36 33

	
37 34
namespace lemon {
38 35

	
... ...
@@ -50,14 +47,19 @@
50 47
  ///
51 48
  /// In general this class is the fastest implementation available
52 49
  /// in LEMON for the minimum cost flow problem.
53
  /// Moreover it supports both direction of the supply/demand inequality
54
  /// constraints. For more information see \ref ProblemType.
50
  /// Moreover it supports both directions of the supply/demand inequality
51
  /// constraints. For more information see \ref SupplyType.
52
  ///
53
  /// Most of the parameters of the problem (except for the digraph)
54
  /// can be given using separate functions, and the algorithm can be
55
  /// executed using the \ref run() function. If some parameters are not
56
  /// specified, then default values will be used.
55 57
  ///
56 58
  /// \tparam GR The digraph type the algorithm runs on.
57
  /// \tparam F The value type used for flow amounts, capacity bounds
59
  /// \tparam V The value type used for flow amounts, capacity bounds
58 60
  /// and supply values in the algorithm. By default it is \c int.
59 61
  /// \tparam C The value type used for costs and potentials in the
60
  /// algorithm. By default it is the same as \c 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
63 65
  /// be integer.
... ...
@@ -65,34 +67,92 @@
65 67
  /// \note %NetworkSimplex provides five different pivot rule
66 68
  /// implementations, from which the most efficient one is used
67 69
  /// by default. For more information see \ref PivotRule.
68
  template <typename GR, typename 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 algorithm
74
    typedef F Flow;
75
    /// The cost type of the algorithm
75
    /// The type of the flow amounts, capacity bounds and supply values
76
    typedef V Value;
77
    /// The type of the arc costs
76 78
    typedef C Cost;
77
#ifdef DOXYGEN
78
    /// The type of the flow map
79
    typedef GR::ArcMap<Flow> FlowMap;
80
    /// The type of the potential map
81
    typedef GR::NodeMap<Cost> PotentialMap;
82
#else
83
    /// The type of the flow map
84
    typedef typename GR::template ArcMap<Flow> FlowMap;
85
    /// The type of the potential map
86
    typedef typename GR::template NodeMap<Cost> PotentialMap;
87
#endif
88 79

	
89 80
  public:
90 81

	
91
    /// \brief Enum type for selecting the pivot rule.
82
    /// \brief Problem type constants for the \c run() function.
92 83
    ///
93
    /// Enum type for selecting the pivot rule for the \ref run()
84
    /// Enum type containing the problem type constants that can be
85
    /// returned by the \ref run() function of the algorithm.
86
    enum ProblemType {
87
      /// The problem has no feasible solution (flow).
88
      INFEASIBLE,
89
      /// The problem has optimal solution (i.e. it is feasible and
90
      /// bounded), and the algorithm has found optimal flow and node
91
      /// potentials (primal and dual solutions).
92
      OPTIMAL,
93
      /// The objective function of the problem is unbounded, i.e.
94
      /// there is a directed cycle having negative total cost and
95
      /// infinite upper bound.
96
      UNBOUNDED
97
    };
98
    
99
    /// \brief Constants for selecting the type of the supply constraints.
100
    ///
101
    /// Enum type containing constants for selecting the supply type,
102
    /// i.e. the direction of the inequalities in the supply/demand
103
    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
104
    ///
105
    /// The default supply type is \c GEQ, since this form is supported
106
    /// by other minimum cost flow algorithms and the \ref Circulation
107
    /// algorithm, as well.
108
    /// The \c LEQ problem type can be selected using the \ref supplyType()
94 109
    /// function.
95 110
    ///
111
    /// Note that the equality form is a special case of both supply types.
112
    enum SupplyType {
113

	
114
      /// This option means that there are <em>"greater or equal"</em>
115
      /// supply/demand constraints in the definition, i.e. the exact
116
      /// formulation of the problem is the following.
117
      /**
118
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
119
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
120
              sup(u) \quad \forall u\in V \f]
121
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
122
      */
123
      /// It means that the total demand must be greater or equal to the 
124
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
125
      /// negative) and all the supplies have to be carried out from 
126
      /// the supply nodes, but there could be demands that are not 
127
      /// satisfied.
128
      GEQ,
129
      /// It is just an alias for the \c GEQ option.
130
      CARRY_SUPPLIES = GEQ,
131

	
132
      /// This option means that there are <em>"less or equal"</em>
133
      /// supply/demand constraints in the definition, i.e. the exact
134
      /// formulation of the problem is the following.
135
      /**
136
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
137
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
138
              sup(u) \quad \forall u\in V \f]
139
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
140
      */
141
      /// It means that the total demand must be less or equal to the 
142
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
143
      /// positive) and all the demands have to be satisfied, but there
144
      /// could be supplies that are not carried out from the supply
145
      /// nodes.
146
      LEQ,
147
      /// It is just an alias for the \c LEQ option.
148
      SATISFY_DEMANDS = LEQ
149
    };
150
    
151
    /// \brief Constants for selecting the pivot rule.
152
    ///
153
    /// Enum type containing constants for selecting the pivot rule for
154
    /// the \ref run() function.
155
    ///
96 156
    /// \ref NetworkSimplex provides five different pivot rule
97 157
    /// implementations that significantly affect the running time
98 158
    /// of the algorithm.
... ...
@@ -131,71 +191,15 @@
131 191
      ALTERING_LIST
132 192
    };
133 193
    
134
    /// \brief Enum type for selecting the problem type.
135
    ///
136
    /// Enum type for selecting the problem type, i.e. the direction of
137
    /// the inequalities in the supply/demand constraints of the
138
    /// \ref min_cost_flow "minimum cost flow problem".
139
    ///
140
    /// The default problem type is \c GEQ, since this form is supported
141
    /// by other minimum cost flow algorithms and the \ref Circulation
142
    /// algorithm as well.
143
    /// The \c LEQ problem type can be selected using the \ref problemType()
144
    /// function.
145
    ///
146
    /// Note that the equality form is a special case of both problem type.
147
    enum ProblemType {
148

	
149
      /// This option means that there are "<em>greater or equal</em>"
150
      /// constraints in the defintion, i.e. the exact formulation of the
151
      /// problem is the following.
152
      /**
153
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
154
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
155
              sup(u) \quad \forall u\in V \f]
156
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
157
      */
158
      /// It means that the total demand must be greater or equal to the 
159
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
160
      /// negative) and all the supplies have to be carried out from 
161
      /// the supply nodes, but there could be demands that are not 
162
      /// satisfied.
163
      GEQ,
164
      /// It is just an alias for the \c GEQ option.
165
      CARRY_SUPPLIES = GEQ,
166

	
167
      /// This option means that there are "<em>less or equal</em>"
168
      /// constraints in the defintion, i.e. the exact formulation of the
169
      /// problem is the following.
170
      /**
171
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
172
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
173
              sup(u) \quad \forall u\in V \f]
174
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
175
      */
176
      /// It means that the total demand must be less or equal to the 
177
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
178
      /// positive) and all the demands have to be satisfied, but there
179
      /// could be supplies that are not carried out from the supply
180
      /// nodes.
181
      LEQ,
182
      /// It is just an alias for the \c LEQ option.
183
      SATISFY_DEMANDS = LEQ
184
    };
185

	
186 194
  private:
187 195

	
188 196
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
189 197

	
190
    typedef typename GR::template ArcMap<Flow> FlowArcMap;
191
    typedef typename GR::template ArcMap<Cost> CostArcMap;
192
    typedef typename GR::template NodeMap<Flow> FlowNodeMap;
193

	
194 198
    typedef std::vector<Arc> ArcVector;
195 199
    typedef std::vector<Node> NodeVector;
196 200
    typedef std::vector<int> IntVector;
197 201
    typedef std::vector<bool> BoolVector;
198
    typedef std::vector<Flow> FlowVector;
202
    typedef std::vector<Value> ValueVector;
199 203
    typedef std::vector<Cost> CostVector;
200 204

	
201 205
    // State constants for arcs
... ...
@@ -213,32 +217,23 @@
213 217
    int _arc_num;
214 218

	
215 219
    // Parameters of the problem
216
    FlowArcMap *_plower;
217
    FlowArcMap *_pupper;
218
    CostArcMap *_pcost;
219
    FlowNodeMap *_psupply;
220
    bool _pstsup;
221
    Node _psource, _ptarget;
222
    Flow _pstflow;
223
    ProblemType _ptype;
224

	
225
    // Result maps
226
    FlowMap *_flow_map;
227
    PotentialMap *_potential_map;
228
    bool _local_flow;
229
    bool _local_potential;
220
    bool _have_lower;
221
    SupplyType _stype;
222
    Value _sum_supply;
230 223

	
231 224
    // Data structures for storing the digraph
232 225
    IntNodeMap _node_id;
233
    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

	
244 239
    // Data for storing the spanning tree structure
... ...
@@ -257,7 +252,16 @@
257 252
    int in_arc, join, u_in, v_in, u_out, v_out;
258 253
    int first, second, right, last;
259 254
    int stem, par_stem, new_stem;
260
    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:
263 267

	
... ...
@@ -659,25 +663,68 @@
659 663
    ///
660 664
    /// \param graph The digraph the algorithm runs on.
661 665
    NetworkSimplex(const GR& graph) :
662
      _graph(graph),
663
      _plower(NULL), _pupper(NULL), _pcost(NULL),
664
      _psupply(NULL), _pstsup(false), _ptype(GEQ),
665
      _flow_map(NULL), _potential_map(NULL),
666
      _local_flow(false), _local_potential(false),
667
      _node_id(graph)
666
      _graph(graph), _node_id(graph), _arc_id(graph),
667
      INF(std::numeric_limits<Value>::has_infinity ?
668
          std::numeric_limits<Value>::infinity() :
669
          std::numeric_limits<Value>::max())
668 670
    {
669
      LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
670
                   std::numeric_limits<Flow>::is_signed,
671
        "The flow type of NetworkSimplex must be signed integer");
672
      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
673
                   std::numeric_limits<Cost>::is_signed,
674
        "The cost type of NetworkSimplex must be signed integer");
671
      // Check the value types
672
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
673
        "The flow type of NetworkSimplex must be signed");
674
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
675
        "The cost type of NetworkSimplex must be signed");
676
        
677
      // Resize vectors
678
      _node_num = countNodes(_graph);
679
      _arc_num = countArcs(_graph);
680
      int all_node_num = _node_num + 1;
681
      int all_arc_num = _arc_num + _node_num;
682

	
683
      _source.resize(all_arc_num);
684
      _target.resize(all_arc_num);
685

	
686
      _lower.resize(all_arc_num);
687
      _upper.resize(all_arc_num);
688
      _cap.resize(all_arc_num);
689
      _cost.resize(all_arc_num);
690
      _supply.resize(all_node_num);
691
      _flow.resize(all_arc_num);
692
      _pi.resize(all_node_num);
693

	
694
      _parent.resize(all_node_num);
695
      _pred.resize(all_node_num);
696
      _forward.resize(all_node_num);
697
      _thread.resize(all_node_num);
698
      _rev_thread.resize(all_node_num);
699
      _succ_num.resize(all_node_num);
700
      _last_succ.resize(all_node_num);
701
      _state.resize(all_arc_num);
702

	
703
      // Copy the graph (store the arcs in a mixed order)
704
      int i = 0;
705
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
706
        _node_id[n] = i;
707
      }
708
      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
709
      i = 0;
710
      for (ArcIt a(_graph); a != INVALID; ++a) {
711
        _arc_id[a] = i;
712
        _source[i] = _node_id[_graph.source(a)];
713
        _target[i] = _node_id[_graph.target(a)];
714
        if ((i += k) >= _arc_num) i = (i % k) + 1;
675 715
    }
676 716

	
677
    /// Destructor.
678
    ~NetworkSimplex() {
679
      if (_local_flow) delete _flow_map;
680
      if (_local_potential) delete _potential_map;
717
      // Initialize maps
718
      for (int i = 0; i != _node_num; ++i) {
719
        _supply[i] = 0;
720
      }
721
      for (int i = 0; i != _arc_num; ++i) {
722
        _lower[i] = 0;
723
        _upper[i] = INF;
724
        _cost[i] = 1;
725
      }
726
      _have_lower = false;
727
      _stype = GEQ;
681 728
    }
682 729

	
683 730
    /// \name Parameters
... ...
@@ -689,21 +736,19 @@
689 736
    /// \brief Set the lower bounds on the arcs.
690 737
    ///
691 738
    /// This function sets the lower bounds on the arcs.
692
    /// If neither this function nor \ref boundMaps() is used before
693
    /// calling \ref run(), the lower bounds will be set to zero
694
    /// 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 Flow type
743
    /// Its \c Value type must be convertible to the \c Value type
698 744
    /// of the algorithm.
699 745
    ///
700 746
    /// \return <tt>(*this)</tt>
701
    template <typename LOWER>
702
    NetworkSimplex& lowerMap(const LOWER& map) {
703
      delete _plower;
704
      _plower = new FlowArcMap(_graph);
747
    template <typename LowerMap>
748
    NetworkSimplex& lowerMap(const LowerMap& map) {
749
      _have_lower = true;
705 750
      for (ArcIt a(_graph); a != INVALID; ++a) {
706
        (*_plower)[a] = map[a];
751
        _lower[_arc_id[a]] = map[a];
707 752
      }
708 753
      return *this;
709 754
    }
... ...
@@ -711,63 +756,23 @@
711 756
    /// \brief Set the upper bounds (capacities) on the arcs.
712 757
    ///
713 758
    /// This function sets the upper bounds (capacities) on the arcs.
714
    /// If none of the functions \ref upperMap(), \ref capacityMap()
715
    /// and \ref boundMaps() is used before calling \ref run(),
716
    /// the upper bounds (capacities) will be set to
717
    /// \c std::numeric_limits<Flow>::max() on all arcs.
759
    /// If it is not used before calling \ref run(), the upper bounds
760
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
761
    /// unbounded from above on each arc).
718 762
    ///
719 763
    /// \param map An arc map storing the upper bounds.
720
    /// Its \c Value type must be convertible to the \c Flow type
764
    /// Its \c Value type must be convertible to the \c Value type
721 765
    /// of the algorithm.
722 766
    ///
723 767
    /// \return <tt>(*this)</tt>
724
    template<typename UPPER>
725
    NetworkSimplex& upperMap(const UPPER& map) {
726
      delete _pupper;
727
      _pupper = new FlowArcMap(_graph);
768
    template<typename UpperMap>
769
    NetworkSimplex& upperMap(const UpperMap& map) {
728 770
      for (ArcIt a(_graph); a != INVALID; ++a) {
729
        (*_pupper)[a] = map[a];
771
        _upper[_arc_id[a]] = map[a];
730 772
      }
731 773
      return *this;
732 774
    }
733 775

	
734
    /// \brief Set the upper bounds (capacities) on the arcs.
735
    ///
736
    /// This function sets the upper bounds (capacities) on the arcs.
737
    /// It is just an alias for \ref upperMap().
738
    ///
739
    /// \return <tt>(*this)</tt>
740
    template<typename CAP>
741
    NetworkSimplex& capacityMap(const CAP& map) {
742
      return upperMap(map);
743
    }
744

	
745
    /// \brief Set the lower and upper bounds on the arcs.
746
    ///
747
    /// This function sets the lower and upper bounds on the arcs.
748
    /// If neither this function nor \ref lowerMap() is used before
749
    /// calling \ref run(), the lower bounds will be set to zero
750
    /// on all arcs.
751
    /// If none of the functions \ref upperMap(), \ref capacityMap()
752
    /// and \ref boundMaps() is used before calling \ref run(),
753
    /// the upper bounds (capacities) will be set to
754
    /// \c std::numeric_limits<Flow>::max() on all arcs.
755
    ///
756
    /// \param lower An arc map storing the lower bounds.
757
    /// \param upper An arc map storing the upper bounds.
758
    ///
759
    /// The \c Value type of the maps must be convertible to the
760
    /// \c Flow type of the algorithm.
761
    ///
762
    /// \note This function is just a shortcut of calling \ref lowerMap()
763
    /// and \ref upperMap() separately.
764
    ///
765
    /// \return <tt>(*this)</tt>
766
    template <typename LOWER, typename UPPER>
767
    NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
768
      return lowerMap(lower).upperMap(upper);
769
    }
770

	
771 776
    /// \brief Set the costs of the arcs.
772 777
    ///
773 778
    /// This function sets the costs of the arcs.
... ...
@@ -779,12 +784,10 @@
779 784
    /// of the algorithm.
780 785
    ///
781 786
    /// \return <tt>(*this)</tt>
782
    template<typename COST>
783
    NetworkSimplex& costMap(const COST& map) {
784
      delete _pcost;
785
      _pcost = new CostArcMap(_graph);
787
    template<typename CostMap>
788
    NetworkSimplex& costMap(const CostMap& map) {
786 789
      for (ArcIt a(_graph); a != INVALID; ++a) {
787
        (*_pcost)[a] = map[a];
790
        _cost[_arc_id[a]] = map[a];
788 791
      }
789 792
      return *this;
790 793
    }
... ...
@@ -797,17 +800,14 @@
797 800
    /// (It makes sense only if non-zero lower bounds are given.)
798 801
    ///
799 802
    /// \param map A node map storing the supply values.
800
    /// Its \c Value type must be convertible to the \c Flow type
803
    /// Its \c Value type must be convertible to the \c Value type
801 804
    /// of the algorithm.
802 805
    ///
803 806
    /// \return <tt>(*this)</tt>
804
    template<typename SUP>
805
    NetworkSimplex& supplyMap(const SUP& map) {
806
      delete _psupply;
807
      _pstsup = false;
808
      _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;
813 813
    }
... ...
@@ -820,68 +820,36 @@
820 820
    /// calling \ref run(), the supply of each node will be set to zero.
821 821
    /// (It makes sense only if non-zero lower bounds are given.)
822 822
    ///
823
    /// Using this function has the same effect as using \ref supplyMap()
824
    /// with such a map in which \c k is assigned to \c s, \c -k is
825
    /// assigned to \c t and all other nodes have zero supply value.
826
    ///
823 827
    /// \param s The source node.
824 828
    /// \param t The target node.
825 829
    /// \param k The required amount of flow from node \c s to node \c t
826 830
    /// (i.e. the supply of \c s and the demand of \c t).
827 831
    ///
828 832
    /// \return <tt>(*this)</tt>
829
    NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
830
      delete _psupply;
831
      _psupply = NULL;
832
      _pstsup = true;
833
      _psource = s;
834
      _ptarget = t;
835
      _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.
842
    /// \brief Set the type of the supply constraints.
840 843
    ///
841
    /// This function sets the problem type for the algorithm.
842
    /// If it is not used before calling \ref run(), the \ref GEQ problem
844
    /// This function sets the type of the supply/demand constraints.
845
    /// If it is not used before calling \ref run(), the \ref GEQ supply
843 846
    /// type will be used.
844 847
    ///
845
    /// For more information see \ref 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;
850
      return *this;
851
    }
852

	
853
    /// \brief Set the flow map.
854
    ///
855
    /// This function sets the flow map.
856
    /// If it is not used before calling \ref run(), an instance will
857
    /// be allocated automatically. The destructor deallocates this
858
    /// automatically allocated map, of course.
859
    ///
860
    /// \return <tt>(*this)</tt>
861
    NetworkSimplex& flowMap(FlowMap& map) {
862
      if (_local_flow) {
863
        delete _flow_map;
864
        _local_flow = false;
865
      }
866
      _flow_map = &map;
867
      return *this;
868
    }
869

	
870
    /// \brief Set the potential map.
871
    ///
872
    /// This function sets the potential map, which is used for storing
873
    /// the dual solution.
874
    /// If it is not used before calling \ref run(), an instance will
875
    /// be allocated automatically. The destructor deallocates this
876
    /// automatically allocated map, of course.
877
    ///
878
    /// \return <tt>(*this)</tt>
879
    NetworkSimplex& potentialMap(PotentialMap& map) {
880
      if (_local_potential) {
881
        delete _potential_map;
882
        _local_potential = false;
883
      }
884
      _potential_map = &map;
851
    NetworkSimplex& supplyType(SupplyType supply_type) {
852
      _stype = supply_type;
885 853
      return *this;
886 854
    }
887 855
    
... ...
@@ -896,13 +864,12 @@
896 864
    ///
897 865
    /// This function runs the algorithm.
898 866
    /// The paramters can be specified using functions \ref lowerMap(),
899
    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
900
    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
901
    /// \ref 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
908 875
    ///
... ...
@@ -910,33 +877,44 @@
910 877
    /// that have been given are kept for the next call, unless
911 878
    /// \ref reset() is called, thus only the modified parameters
912 879
    /// have to be set again. See \ref reset() for examples.
880
    /// However the underlying digraph must not be modified after this
881
    /// class have been constructed, since it copies and extends the graph.
913 882
    ///
914 883
    /// \param pivot_rule The pivot rule that will be used during the
915 884
    /// algorithm. For more information see \ref PivotRule.
916 885
    ///
917
    /// \return \c true if a feasible flow can be found.
918
    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
919
      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

	
922 900
    /// \brief Reset all the parameters that have been given before.
923 901
    ///
924 902
    /// This function resets all the paramaters that have been given
925 903
    /// before using functions \ref lowerMap(), \ref upperMap(),
926
    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
927
    /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
928
    /// \ref 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,
935 913
    /// \code
936 914
    ///   NetworkSimplex<ListDigraph> ns(graph);
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
    ///
942 920
    ///   // Run again with modified cost map (reset() is not called,
... ...
@@ -947,29 +925,22 @@
947 925
    ///   // Run again from scratch using reset()
948 926
    ///   // (the lower bounds will be set to zero on all arcs)
949 927
    ///   ns.reset();
950
    ///   ns.capacityMap(cap).costMap(cost)
928
    ///   ns.upperMap(capacity).costMap(cost)
951 929
    ///     .supplyMap(sup).run();
952 930
    /// \endcode
953 931
    ///
954 932
    /// \return <tt>(*this)</tt>
955 933
    NetworkSimplex& reset() {
956
      delete _plower;
957
      delete _pupper;
958
      delete _pcost;
959
      delete _psupply;
960
      _plower = NULL;
961
      _pupper = NULL;
962
      _pcost = NULL;
963
      _psupply = NULL;
964
      _pstsup = false;
965
      _ptype = GEQ;
966
      if (_local_flow) delete _flow_map;
967
      if (_local_potential) delete _potential_map;
968
      _flow_map = NULL;
969
      _potential_map = NULL;
970
      _local_flow = false;
971
      _local_potential = false;
972

	
934
      for (int i = 0; i != _node_num; ++i) {
935
        _supply[i] = 0;
936
      }
937
      for (int i = 0; i != _arc_num; ++i) {
938
        _lower[i] = 0;
939
        _upper[i] = INF;
940
        _cost[i] = 1;
941
      }
942
      _have_lower = false;
943
      _stype = GEQ;
973 944
      return *this;
974 945
    }
975 946

	
... ...
@@ -985,7 +956,7 @@
985 956
    /// \brief Return the total cost of the found flow.
986 957
    ///
987 958
    /// This function returns the total cost of the found flow.
988
    /// The complexity of the function is O(e).
959
    /// Its complexity is O(e).
989 960
    ///
990 961
    /// \note The return type of the function can be specified as a
991 962
    /// template parameter. For example,
... ...
@@ -997,15 +968,12 @@
997 968
    /// function.
998 969
    ///
999 970
    /// \pre \ref run() must be called before using this function.
1000
    template <typename Num>
1001
    Num totalCost() const {
1002
      Num c = 0;
1003
      if (_pcost) {
1004
        for (ArcIt e(_graph); e != INVALID; ++e)
1005
          c += (*_flow_map)[e] * (*_pcost)[e];
1006
      } else {
1007
        for (ArcIt e(_graph); e != INVALID; ++e)
1008
          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;
1011 979
    }
... ...
@@ -1021,18 +989,22 @@
1021 989
    /// This function returns the flow on the given arc.
1022 990
    ///
1023 991
    /// \pre \ref run() must be called before using this function.
1024
    Flow flow(const Arc& a) const {
1025
      return (*_flow_map)[a];
992
    Value flow(const Arc& a) const {
993
      return _flow[_arc_id[a]];
1026 994
    }
1027 995

	
1028
    /// \brief Return a const reference to the flow map.
996
    /// \brief Return the flow map (the primal solution).
1029 997
    ///
1030
    /// This function returns a const reference to an arc map storing
1031
    /// the found flow.
998
    /// This function copies the flow value on each arc into the given
999
    /// map. The \c Value type of the algorithm must be convertible to
1000
    /// the \c Value type of the map.
1032 1001
    ///
1033 1002
    /// \pre \ref run() must be called before using this function.
1034
    const FlowMap& flowMap() const {
1035
      return *_flow_map;
1003
    template <typename FlowMap>
1004
    void flowMap(FlowMap &map) const {
1005
      for (ArcIt a(_graph); a != INVALID; ++a) {
1006
        map.set(a, _flow[_arc_id[a]]);
1007
      }
1036 1008
    }
1037 1009

	
1038 1010
    /// \brief Return the potential (dual value) of the given node.
... ...
@@ -1042,19 +1014,22 @@
1042 1014
    ///
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];
1017
      return _pi[_node_id[n]];
1046 1018
    }
1047 1019

	
1048
    /// \brief Return a const reference to the potential map
1049
    /// (the dual solution).
1020
    /// \brief Return the potential map (the dual solution).
1050 1021
    ///
1051
    /// This function returns a const reference to a node map storing
1052
    /// the found potentials, which form the dual solution of the
1053
    /// \ref min_cost_flow "minimum cost flow" problem.
1022
    /// This function copies the potential (dual value) of each node
1023
    /// into the given map.
1024
    /// The \c Cost type of the algorithm must be convertible to the
1025
    /// \c Value type of the map.
1054 1026
    ///
1055 1027
    /// \pre \ref run() must be called before using this function.
1056
    const PotentialMap& potentialMap() const {
1057
      return *_potential_map;
1028
    template <typename PotentialMap>
1029
    void potentialMap(PotentialMap &map) const {
1030
      for (NodeIt n(_graph); n != INVALID; ++n) {
1031
        map.set(n, _pi[_node_id[n]]);
1032
      }
1058 1033
    }
1059 1034

	
1060 1035
    /// @}
... ...
@@ -1063,155 +1038,51 @@
1063 1038

	
1064 1039
    // Initialize internal data structures
1065 1040
    bool init() {
1066
      // Initialize result maps
1067
      if (!_flow_map) {
1068
        _flow_map = new FlowMap(_graph);
1069
        _local_flow = true;
1041
      if (_node_num == 0) return false;
1042

	
1043
      // Check the sum of supply values
1044
      _sum_supply = 0;
1045
      for (int i = 0; i != _node_num; ++i) {
1046
        _sum_supply += _supply[i];
1070 1047
      }
1071
      if (!_potential_map) {
1072
        _potential_map = new PotentialMap(_graph);
1073
        _local_potential = true;
1048
      if ( !((_stype == GEQ && _sum_supply <= 0) ||
1049
             (_stype == LEQ && _sum_supply >= 0)) ) return false;
1050

	
1051
      // Remove non-zero lower bounds
1052
      if (_have_lower) {
1053
        for (int i = 0; i != _arc_num; ++i) {
1054
          Value c = _lower[i];
1055
          if (c >= 0) {
1056
            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
1057
          } else {
1058
            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
1059
          }
1060
          _supply[_source[i]] -= c;
1061
          _supply[_target[i]] += c;
1062
        }
1063
      } else {
1064
        for (int i = 0; i != _arc_num; ++i) {
1065
          _cap[i] = _upper[i];
1066
        }
1074 1067
      }
1075 1068

	
1076
      // Initialize vectors
1077
      _node_num = countNodes(_graph);
1078
      _arc_num = countArcs(_graph);
1079
      int all_node_num = _node_num + 1;
1080
      int all_arc_num = _arc_num + _node_num;
1081
      if (_node_num == 0) return false;
1082

	
1083
      _arc_ref.resize(_arc_num);
1084
      _source.resize(all_arc_num);
1085
      _target.resize(all_arc_num);
1086

	
1087
      _cap.resize(all_arc_num);
1088
      _cost.resize(all_arc_num);
1089
      _supply.resize(all_node_num);
1090
      _flow.resize(all_arc_num);
1091
      _pi.resize(all_node_num);
1092

	
1093
      _parent.resize(all_node_num);
1094
      _pred.resize(all_node_num);
1095
      _forward.resize(all_node_num);
1096
      _thread.resize(all_node_num);
1097
      _rev_thread.resize(all_node_num);
1098
      _succ_num.resize(all_node_num);
1099
      _last_succ.resize(all_node_num);
1100
      _state.resize(all_arc_num);
1101

	
1102
      // Initialize node related data
1103
      bool valid_supply = true;
1104
      Flow sum_supply = 0;
1105
      if (!_pstsup && !_psupply) {
1106
        _pstsup = true;
1107
        _psource = _ptarget = NodeIt(_graph);
1108
        _pstflow = 0;
1069
      // Initialize artifical cost
1070
      Cost ART_COST;
1071
      if (std::numeric_limits<Cost>::is_exact) {
1072
        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
1073
      } else {
1074
        ART_COST = std::numeric_limits<Cost>::min();
1075
        for (int i = 0; i != _arc_num; ++i) {
1076
          if (_cost[i] > ART_COST) ART_COST = _cost[i];
1109 1077
      }
1110
      if (_psupply) {
1111
        int i = 0;
1112
        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1113
          _node_id[n] = i;
1114
          _supply[i] = (*_psupply)[n];
1115
          sum_supply += _supply[i];
1116
        }
1117
        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
1118
                       (_ptype == LEQ && sum_supply >= 0);
1119
      } else {
1120
        int i = 0;
1121
        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1122
          _node_id[n] = i;
1123
          _supply[i] = 0;
1124
        }
1125
        _supply[_node_id[_psource]] =  _pstflow;
1126
        _supply[_node_id[_ptarget]] = -_pstflow;
1127
      }
1128
      if (!valid_supply) return false;
1129

	
1130
      // Infinite capacity value
1131
      Flow inf_cap =
1132
        std::numeric_limits<Flow>::has_infinity ?
1133
        std::numeric_limits<Flow>::infinity() :
1134
        std::numeric_limits<Flow>::max();
1135

	
1136
      // Initialize artifical cost
1137
      Cost art_cost;
1138
      if (std::numeric_limits<Cost>::is_exact) {
1139
        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
1140
      } else {
1141
        art_cost = std::numeric_limits<Cost>::min();
1142
        for (int i = 0; i != _arc_num; ++i) {
1143
          if (_cost[i] > art_cost) art_cost = _cost[i];
1144
        }
1145
        art_cost = (art_cost + 1) * _node_num;
1078
        ART_COST = (ART_COST + 1) * _node_num;
1146 1079
      }
1147 1080

	
1148
      // Run Circulation to check if a feasible solution exists
1149
      typedef ConstMap<Arc, Flow> ConstArcMap;
1150
      ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
1151
      FlowNodeMap *csup = NULL;
1152
      bool local_csup = false;
1153
      if (_psupply) {
1154
        csup = _psupply;
1155
      } else {
1156
        csup = new FlowNodeMap(_graph, 0);
1157
        (*csup)[_psource] =  _pstflow;
1158
        (*csup)[_ptarget] = -_pstflow;
1159
        local_csup = true;
1081
      // Initialize arc maps
1082
      for (int i = 0; i != _arc_num; ++i) {
1083
        _flow[i] = 0;
1084
        _state[i] = STATE_LOWER;
1160 1085
      }
1161
      bool circ_result = false;
1162
      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
1163
        // GEQ problem type
1164
        if (_plower) {
1165
          if (_pupper) {
1166
            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
1167
              circ(_graph, *_plower, *_pupper, *csup);
1168
            circ_result = circ.run();
1169
          } else {
1170
            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
1171
              circ(_graph, *_plower, inf_arc_map, *csup);
1172
            circ_result = circ.run();
1173
          }
1174
        } else {
1175
          if (_pupper) {
1176
            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
1177
              circ(_graph, zero_arc_map, *_pupper, *csup);
1178
            circ_result = circ.run();
1179
          } else {
1180
            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
1181
              circ(_graph, zero_arc_map, inf_arc_map, *csup);
1182
            circ_result = circ.run();
1183
          }
1184
        }
1185
      } else {
1186
        // LEQ problem type
1187
        typedef ReverseDigraph<const GR> RevGraph;
1188
        typedef NegMap<FlowNodeMap> NegNodeMap;
1189
        RevGraph rgraph(_graph);
1190
        NegNodeMap neg_csup(*csup);
1191
        if (_plower) {
1192
          if (_pupper) {
1193
            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
1194
              circ(rgraph, *_plower, *_pupper, neg_csup);
1195
            circ_result = circ.run();
1196
          } else {
1197
            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
1198
              circ(rgraph, *_plower, inf_arc_map, neg_csup);
1199
            circ_result = circ.run();
1200
          }
1201
        } else {
1202
          if (_pupper) {
1203
            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
1204
              circ(rgraph, zero_arc_map, *_pupper, neg_csup);
1205
            circ_result = circ.run();
1206
          } else {
1207
            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
1208
              circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
1209
            circ_result = circ.run();
1210
          }
1211
        }
1212
      }
1213
      if (local_csup) delete csup;
1214
      if (!circ_result) return false;
1215 1086

	
1216 1087
      // Set data for the artificial root node
1217 1088
      _root = _node_num;
... ...
@@ -1219,89 +1090,30 @@
1219 1090
      _pred[_root] = -1;
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
      }
1307 1119

	
... ...
@@ -1336,13 +1148,14 @@
1336 1148
      }
1337 1149
      delta = _cap[in_arc];
1338 1150
      int result = 0;
1339
      Flow d;
1151
      Value d;
1340 1152
      int e;
1341 1153

	
1342 1154
      // Search the cycle along the path form the first node to the root
1343 1155
      for (int u = first; u != join; u = _parent[u]) {
1344 1156
        e = _pred[u];
1345
        d = _forward[u] ? _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;
1348 1161
          u_out = u;
... ...
@@ -1352,7 +1165,8 @@
1352 1165
      // Search the cycle along the path form the second node to the root
1353 1166
      for (int u = second; u != join; u = _parent[u]) {
1354 1167
        e = _pred[u];
1355
        d = _forward[u] ? _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;
1358 1172
          u_out = u;
... ...
@@ -1374,7 +1188,7 @@
1374 1188
    void changeFlow(bool change) {
1375 1189
      // Augment along the cycle
1376 1190
      if (delta > 0) {
1377
        Flow val = _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]) {
1380 1194
          _flow[_pred[u]] += _forward[u] ? -val : val;
... ...
@@ -1526,7 +1340,7 @@
1526 1340
    }
1527 1341

	
1528 1342
    // Execute the algorithm
1529
    bool start(PivotRule pivot_rule) {
1343
    ProblemType start(PivotRule pivot_rule) {
1530 1344
      // Select the pivot rule implementation
1531 1345
      switch (pivot_rule) {
1532 1346
        case FIRST_ELIGIBLE:
... ...
@@ -1540,17 +1354,18 @@
1540 1354
        case ALTERING_LIST:
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
    bool start() {
1361
    ProblemType start() {
1548 1362
      PivotRuleImpl pivot(*this);
1549 1363

	
1550 1364
      // Execute the Network Simplex algorithm
1551 1365
      while (pivot.findEnteringArc()) {
1552 1366
        findJoinNode();
1553 1367
        bool change = findLeavingArc();
1368
        if (delta >= INF) return UNBOUNDED;
1554 1369
        changeFlow(change);
1555 1370
        if (change) {
1556 1371
          updateTreeStructure();
... ...
@@ -1558,23 +1373,36 @@
1558 1373
        }
1559 1374
      }
1560 1375

	
1561
      // Copy flow values to _flow_map
1562
      if (_plower) {
1563
        for (int i = 0; i != _arc_num; ++i) {
1564
          Arc e = _arc_ref[i];
1565
          _flow_map->set(e, (*_plower)[e] + _flow[i]);
1376
      // Check feasibility
1377
      if (_sum_supply < 0) {
1378
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1379
          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
1566 1380
        }
1567
      } else {
1568
        for (int i = 0; i != _arc_num; ++i) {
1569
          _flow_map->set(_arc_ref[i], _flow[i]);
1570 1381
        }
1382
      else if (_sum_supply > 0) {
1383
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1384
          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
1571 1385
      }
1572
      // Copy potential values to _potential_map
1573
      for (NodeIt n(_graph); n != INVALID; ++n) {
1574
        _potential_map->set(n, _pi[_node_id[n]]);
1386
      }
1387
      else {
1388
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1389
          if (_flow[e] != 0) return INFEASIBLE;
1390
        }
1575 1391
      }
1576 1392

	
1577
      return true;
1393
      // Transform the solution and the supply map to the original form
1394
      if (_have_lower) {
1395
        for (int i = 0; i != _arc_num; ++i) {
1396
          Value c = _lower[i];
1397
          if (c != 0) {
1398
            _flow[i] += c;
1399
            _supply[_source[i]] += c;
1400
            _supply[_target[i]] -= c;
1401
          }
1402
        }
1403
      }
1404

	
1405
      return OPTIMAL;
1578 1406
    }
1579 1407

	
1580 1408
  }; //class NetworkSimplex
Show white space 6 line context
... ...
@@ -46,13 +46,13 @@
46 46
    typedef CAP CapacityMap;
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.
52 52
    ///
53 53
    /// The type of the map that stores the flow values.
54 54
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
55
    typedef typename Digraph::template ArcMap<Flow> FlowMap;
55
    typedef typename Digraph::template ArcMap<Value> FlowMap;
56 56

	
57 57
    /// \brief Instantiates a FlowMap.
58 58
    ///
... ...
@@ -84,7 +84,7 @@
84 84
    /// \brief The tolerance used by the algorithm
85 85
    ///
86 86
    /// The tolerance used by the algorithm to handle inexact computation.
87
    typedef lemon::Tolerance<Flow> Tolerance;
87
    typedef lemon::Tolerance<Value> Tolerance;
88 88

	
89 89
  };
90 90

	
... ...
@@ -125,7 +125,7 @@
125 125
    ///The type of the capacity map.
126 126
    typedef typename Traits::CapacityMap CapacityMap;
127 127
    ///The type of the flow values.
128
    typedef typename Traits::Flow Flow;
128
    typedef typename Traits::Value Value;
129 129

	
130 130
    ///The type of the flow map.
131 131
    typedef typename Traits::FlowMap FlowMap;
... ...
@@ -151,7 +151,7 @@
151 151
    Elevator* _level;
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

	
157 157
    Tolerance _tolerance;
... ...
@@ -470,7 +470,7 @@
470 470
      }
471 471

	
472 472
      for (NodeIt n(_graph); n != INVALID; ++n) {
473
        Flow excess = 0;
473
        Value excess = 0;
474 474
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
475 475
          excess += (*_flow)[e];
476 476
        }
... ...
@@ -519,7 +519,7 @@
519 519
      _level->initFinish();
520 520

	
521 521
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
522
        Flow rem = (*_capacity)[e] - (*_flow)[e];
522
        Value rem = (*_capacity)[e] - (*_flow)[e];
523 523
        if (_tolerance.positive(rem)) {
524 524
          Node u = _graph.target(e);
525 525
          if ((*_level)[u] == _level->maxLevel()) continue;
... ...
@@ -531,7 +531,7 @@
531 531
        }
532 532
      }
533 533
      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
534
        Flow rem = (*_flow)[e];
534
        Value rem = (*_flow)[e];
535 535
        if (_tolerance.positive(rem)) {
536 536
          Node v = _graph.source(e);
537 537
          if ((*_level)[v] == _level->maxLevel()) continue;
... ...
@@ -564,11 +564,11 @@
564 564
        int num = _node_num;
565 565

	
566 566
        while (num > 0 && n != INVALID) {
567
          Flow excess = (*_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
            Flow rem = (*_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);
574 574
            if ((*_level)[v] < level) {
... ...
@@ -591,7 +591,7 @@
591 591
          }
592 592

	
593 593
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
594
            Flow rem = (*_flow)[e];
594
            Value rem = (*_flow)[e];
595 595
            if (!_tolerance.positive(rem)) continue;
596 596
            Node v = _graph.source(e);
597 597
            if ((*_level)[v] < level) {
... ...
@@ -637,11 +637,11 @@
637 637

	
638 638
        num = _node_num * 20;
639 639
        while (num > 0 && n != INVALID) {
640
          Flow excess = (*_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
            Flow rem = (*_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);
647 647
            if ((*_level)[v] < level) {
... ...
@@ -664,7 +664,7 @@
664 664
          }
665 665

	
666 666
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
667
            Flow rem = (*_flow)[e];
667
            Value rem = (*_flow)[e];
668 668
            if (!_tolerance.positive(rem)) continue;
669 669
            Node v = _graph.source(e);
670 670
            if ((*_level)[v] < level) {
... ...
@@ -778,12 +778,12 @@
778 778

	
779 779
      Node n;
780 780
      while ((n = _level->highestActive()) != INVALID) {
781
        Flow excess = (*_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
          Flow rem = (*_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);
789 789
          if ((*_level)[v] < level) {
... ...
@@ -806,7 +806,7 @@
806 806
        }
807 807

	
808 808
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
809
          Flow rem = (*_flow)[e];
809
          Value rem = (*_flow)[e];
810 810
          if (!_tolerance.positive(rem)) continue;
811 811
          Node v = _graph.source(e);
812 812
          if ((*_level)[v] < level) {
... ...
@@ -897,18 +897,18 @@
897 897
    ///
898 898
    /// \pre Either \ref run() or \ref init() must be called before
899 899
    /// using this function.
900
    Flow flowValue() const {
900
    Value flowValue() const {
901 901
      return (*_excess)[_target];
902 902
    }
903 903

	
904
    /// \brief Returns the flow on the given arc.
904
    /// \brief Returns the flow value on the given arc.
905 905
    ///
906
    /// Returns the flow on the given arc. This method can
906
    /// Returns the flow value on the given arc. This method can
907 907
    /// be called after the second phase of the algorithm.
908 908
    ///
909 909
    /// \pre Either \ref run() or \ref init() must be called before
910 910
    /// using this function.
911
    Flow flow(const Arc& arc) const {
911
    Value flow(const Arc& arc) const {
912 912
      return (*_flow)[arc];
913 913
    }
914 914

	
Show white space 6 line context
... ...
@@ -21,9 +21,7 @@
21 21
#include "test_tools.h"
22 22
#include <lemon/tolerance.h>
23 23

	
24
#ifdef HAVE_CONFIG_H
25 24
#include <lemon/config.h>
26
#endif
27 25

	
28 26
#ifdef LEMON_HAVE_GLPK
29 27
#include <lemon/glpk.h>
Show white space 6 line context
... ...
@@ -18,6 +18,7 @@
18 18

	
19 19
#include <iostream>
20 20
#include <fstream>
21
#include <limits>
21 22

	
22 23
#include <lemon/list_graph.h>
23 24
#include <lemon/lgf_reader.h>
... ...
@@ -33,57 +34,57 @@
33 34

	
34 35
char test_lgf[] =
35 36
  "@nodes\n"
36
  "label  sup1 sup2 sup3 sup4 sup5\n"
37
  "    1    20   27    0   20   30\n"
38
  "    2    -4    0    0   -8   -3\n"
39
  "    3     0    0    0    0    0\n"
40
  "    4     0    0    0    0    0\n"
41
  "    5     9    0    0    6   11\n"
42
  "    6    -6    0    0   -5   -6\n"
43
  "    7     0    0    0    0    0\n"
44
  "    8     0    0    0    0    3\n"
45
  "    9     3    0    0    0    0\n"
46
  "   10    -2    0    0   -7   -2\n"
47
  "   11     0    0    0  -10    0\n"
48
  "   12   -20  -27    0  -30  -20\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"
49 50
  "\n"
50 51
  "@arcs\n"
51
  "       cost  cap low1 low2\n"
52
  " 1  2    70   11    0    8\n"
53
  " 1  3   150    3    0    1\n"
54
  " 1  4    80   15    0    2\n"
55
  " 2  8    80   12    0    0\n"
56
  " 3  5   140    5    0    3\n"
57
  " 4  6    60   10    0    1\n"
58
  " 4  7    80    2    0    0\n"
59
  " 4  8   110    3    0    0\n"
60
  " 5  7    60   14    0    0\n"
61
  " 5 11   120   12    0    0\n"
62
  " 6  3     0    3    0    0\n"
63
  " 6  9   140    4    0    0\n"
64
  " 6 10    90    8    0    0\n"
65
  " 7  1    30    5    0    0\n"
66
  " 8 12    60   16    0    4\n"
67
  " 9 12    50    6    0    0\n"
68
  "10 12    70   13    0    5\n"
69
  "10  2   100    7    0    0\n"
70
  "10  7    60   10    0    0\n"
71
  "11 10    20   14    0    6\n"
72
  "12 11    30   10    0    0\n"
52
  "       cost  cap low1 low2 low3\n"
53
  " 1  2    70   11    0    8    8\n"
54
  " 1  3   150    3    0    1    0\n"
55
  " 1  4    80   15    0    2    2\n"
56
  " 2  8    80   12    0    0    0\n"
57
  " 3  5   140    5    0    3    1\n"
58
  " 4  6    60   10    0    1    0\n"
59
  " 4  7    80    2    0    0    0\n"
60
  " 4  8   110    3    0    0    0\n"
61
  " 5  7    60   14    0    0    0\n"
62
  " 5 11   120   12    0    0    0\n"
63
  " 6  3     0    3    0    0    0\n"
64
  " 6  9   140    4    0    0    0\n"
65
  " 6 10    90    8    0    0    0\n"
66
  " 7  1    30    5    0    0   -5\n"
67
  " 8 12    60   16    0    4    3\n"
68
  " 9 12    50    6    0    0    0\n"
69
  "10 12    70   13    0    5    2\n"
70
  "10  2   100    7    0    0    0\n"
71
  "10  7    60   10    0    0   -3\n"
72
  "11 10    20   14    0    6  -20\n"
73
  "12 11    30   10    0    0  -10\n"
73 74
  "\n"
74 75
  "@attributes\n"
75 76
  "source 1\n"
76 77
  "target 12\n";
77 78

	
78 79

	
79
enum ProblemType {
80
enum SupplyType {
80 81
  EQ,
81 82
  GEQ,
82 83
  LEQ
83 84
};
84 85

	
85 86
// Check the interface of an MCF algorithm
86
template <typename GR, typename Flow, typename Cost>
87
template <typename GR, typename Value, typename Cost>
87 88
class McfClassConcept
88 89
{
89 90
public:
... ...
@@ -94,53 +95,46 @@
94 95
      checkConcept<concepts::Digraph, GR>();
95 96

	
96 97
      MCF mcf(g);
98
      const MCF& const_mcf = mcf;
97 99

	
98 100
      b = mcf.reset()
99 101
             .lowerMap(lower)
100 102
             .upperMap(upper)
101
             .capacityMap(upper)
102
             .boundMaps(lower, upper)
103 103
             .costMap(cost)
104 104
             .supplyMap(sup)
105 105
             .stSupply(n, n, k)
106
             .flowMap(flow)
107
             .potentialMap(pot)
108 106
             .run();
109 107
      
110
      const MCF& const_mcf = mcf;
111

	
112
      const typename MCF::FlowMap &fm = const_mcf.flowMap();
113
      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
114

	
115
      v = const_mcf.totalCost();
116
      double x = const_mcf.template totalCost<double>();
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 &sup;
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

	
146 140
};
... ...
@@ -151,7 +145,7 @@
151 145
           typename SM, typename FM >
152 146
bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
153 147
                const SM& supply, const FM& flow,
154
                ProblemType type = EQ )
148
                SupplyType type = EQ )
155 149
{
156 150
  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
157 151

	
... ...
@@ -208,21 +202,25 @@
208 202
// Run a minimum cost flow algorithm and check the results
209 203
template < typename MCF, typename GR,
210 204
           typename LM, typename UM,
211
           typename CM, typename SM >
212
void checkMcf( const MCF& mcf, bool mcf_result,
205
           typename CM, typename SM,
206
           typename PT >
207
void checkMcf( const MCF& mcf, PT mcf_result,
213 208
               const GR& gr, const LM& lower, const UM& upper,
214 209
               const CM& cost, const SM& supply,
215
               bool 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
  }
228 226
}
... ...
@@ -231,11 +229,13 @@
231 229
{
232 230
  // Check the interfaces
233 231
  {
234
    typedef int Flow;
235
    typedef int Cost;
236 232
    typedef concepts::Digraph GR;
237
    checkConcept< McfClassConcept<GR, Flow, Cost>,
238
                  NetworkSimplex<GR, Flow, Cost> >();
233
    checkConcept< McfClassConcept<GR, int, int>,
234
                  NetworkSimplex<GR> >();
235
    checkConcept< McfClassConcept<GR, double, double>,
236
                  NetworkSimplex<GR, double> >();
237
    checkConcept< McfClassConcept<GR, int, double>,
238
                  NetworkSimplex<GR, int, double> >();
239 239
  }
240 240

	
241 241
  // Run various MCF tests
... ...
@@ -244,8 +244,8 @@
244 244

	
245 245
  // Read the test digraph
246 246
  Digraph gr;
247
  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
248
  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
247
  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
248
  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
249 249
  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
250 250
  Node v, w;
251 251

	
... ...
@@ -255,15 +255,57 @@
255 255
    .arcMap("cap", u)
256 256
    .arcMap("low1", l1)
257 257
    .arcMap("low2", l2)
258
    .arcMap("low3", l3)
258 259
    .nodeMap("sup1", s1)
259 260
    .nodeMap("sup2", s2)
260 261
    .nodeMap("sup3", s3)
261 262
    .nodeMap("sup4", s4)
262 263
    .nodeMap("sup5", s5)
264
    .nodeMap("sup6", s6)
263 265
    .node("source", v)
264 266
    .node("target", w)
265 267
    .run();
266 268

	
269
  // Build a test digraph for testing negative costs
270
  Digraph ngr;
271
  Node n1 = ngr.addNode();
272
  Node n2 = ngr.addNode();
273
  Node n3 = ngr.addNode();
274
  Node n4 = ngr.addNode();
275
  Node n5 = ngr.addNode();
276
  Node n6 = ngr.addNode();
277
  Node n7 = ngr.addNode();
278
  
279
  Arc a1 = ngr.addArc(n1, n2);
280
  Arc a2 = ngr.addArc(n1, n3);
281
  Arc a3 = ngr.addArc(n2, n4);
282
  Arc a4 = ngr.addArc(n3, n4);
283
  Arc a5 = ngr.addArc(n3, n2);
284
  Arc a6 = ngr.addArc(n5, n3);
285
  Arc a7 = ngr.addArc(n5, n6);
286
  Arc a8 = ngr.addArc(n6, n7);
287
  Arc a9 = ngr.addArc(n7, n5);
288
  
289
  Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
290
  ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
291
  Digraph::NodeMap<int> ns(ngr, 0);
292
  
293
  nl2[a7] =  1000;
294
  nl2[a8] = -1000;
295
  
296
  ns[n1] =  100;
297
  ns[n4] = -100;
298
  
299
  nc[a1] =  100;
300
  nc[a2] =   30;
301
  nc[a3] =   20;
302
  nc[a4] =   80;
303
  nc[a5] =   50;
304
  nc[a6] =   10;
305
  nc[a7] =   80;
306
  nc[a8] =   30;
307
  nc[a9] = -120;
308

	
267 309
  // A. Test NetworkSimplex with the default pivot rule
268 310
  {
269 311
    NetworkSimplex<Digraph> mcf(gr);
... ...
@@ -271,63 +313,77 @@
271 313
    // Check the equality form
272 314
    mcf.upperMap(u).costMap(c);
273 315
    checkMcf(mcf, mcf.supplyMap(s1).run(),
274
             gr, l1, u, c, s1, true,  5240, "#A1");
316
             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
275 317
    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
276
             gr, l1, u, c, s2, true,  7620, "#A2");
318
             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
277 319
    mcf.lowerMap(l2);
278 320
    checkMcf(mcf, mcf.supplyMap(s1).run(),
279
             gr, l2, u, c, s1, true,  5970, "#A3");
321
             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
280 322
    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
281
             gr, l2, u, c, s2, true,  8010, "#A4");
323
             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
282 324
    mcf.reset();
283 325
    checkMcf(mcf, mcf.supplyMap(s1).run(),
284
             gr, l1, cu, cc, s1, true,  74, "#A5");
326
             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
285 327
    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
286
             gr, l2, cu, cc, s2, true,  94, "#A6");
328
             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
287 329
    mcf.reset();
288 330
    checkMcf(mcf, mcf.run(),
289
             gr, l1, cu, cc, s3, true,   0, "#A7");
290
    checkMcf(mcf, mcf.boundMaps(l2, u).run(),
291
             gr, 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(s4);
339
    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
295 340
    checkMcf(mcf, mcf.run(),
296
             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
297
    mcf.problemType(mcf.GEQ);
341
             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
342
    mcf.supplyType(mcf.GEQ);
298 343
    checkMcf(mcf, mcf.lowerMap(l2).run(),
299
             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
300
    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
344
             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
345
    mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
301 346
    checkMcf(mcf, mcf.run(),
302
             gr, l2, u, c, s5, false,    0, "#A11", GEQ);
347
             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
303 348

	
304 349
    // Check the LEQ form
305
    mcf.reset().problemType(mcf.LEQ);
306
    mcf.upperMap(u).costMap(c).supplyMap(s5);
350
    mcf.reset().supplyType(mcf.LEQ);
351
    mcf.upperMap(u).costMap(c).supplyMap(s6);
307 352
    checkMcf(mcf, mcf.run(),
308
             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
353
             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
309 354
    checkMcf(mcf, mcf.lowerMap(l2).run(),
310
             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
311
    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
355
             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
356
    mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
312 357
    checkMcf(mcf, mcf.run(),
313
             gr, l2, u, c, s4, false,    0, "#A14", LEQ);
358
             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
359

	
360
    // Check negative costs
361
    NetworkSimplex<Digraph> nmcf(ngr);
362
    nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
363
    checkMcf(nmcf, nmcf.run(),
364
      ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
365
    checkMcf(nmcf, nmcf.upperMap(nu2).run(),
366
      ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
367
    nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
368
    checkMcf(nmcf, nmcf.run(),
369
      ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
314 370
  }
315 371

	
316 372
  // B. Test NetworkSimplex with each pivot rule
317 373
  {
318 374
    NetworkSimplex<Digraph> mcf(gr);
319
    mcf.supplyMap(s1).costMap(c).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

	
333 389
  return 0;
Show white space 6 line context
... ...
@@ -18,9 +18,7 @@
18 18

	
19 19
#include "test_tools.h"
20 20

	
21
#ifdef HAVE_CONFIG_H
22 21
#include <lemon/config.h>
23
#endif
24 22

	
25 23
#ifdef LEMON_HAVE_CPLEX
26 24
#include <lemon/cplex.h>
Show white space 6 line context
... ...
@@ -119,8 +119,8 @@
119 119

	
120 120
  ti.restart();
121 121
  NetworkSimplex<Digraph, Value> ns(g);
122
  ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
123
  if (sum_sup > 0) ns.problemType(ns.LEQ);
122
  ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup);
123
  if (sum_sup > 0) ns.supplyType(ns.LEQ);
124 124
  if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';
125 125
  ti.restart();
126 126
  bool res = ns.run();
0 comments (0 inline)