↑ Collapse diff ↑
Ignore white space 6 line context
... ...
@@ -19,4 +19,2 @@
19 19

	
20
ADD_DEFINITIONS(-DHAVE_CONFIG_H)
21

	
22 20
IF(MSVC)
... ...
@@ -30,4 +28,2 @@
30 28

	
31
ADD_DEFINITIONS(-DHAVE_CONFIG_H)
32

	
33 29
INCLUDE(CheckTypeSize)
Ignore white space 6 line context
... ...
@@ -13,7 +13,8 @@
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 \
Ignore white space 6 line context
... ...
@@ -3,24 +3,44 @@
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

	
Ignore 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

	
... ...
@@ -13,3 +20,4 @@
13 20
  cplex91.dll
14
  PATHS "C:/ILOG/CPLEX91/bin/x86_win32")
21
  PATHS "C:/ILOG/CPLEX91/bin/x86_win32"
22
)
15 23

	
... ...
@@ -18,2 +26,5 @@
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)
Ignore 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]")
... ...
@@ -5,10 +7,42 @@
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

	
Ignore white space 6 line context
... ...
@@ -354,8 +354,8 @@
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.
... ...
@@ -364,3 +364,3 @@
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.
... ...
@@ -406,3 +406,3 @@
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$
... ...
@@ -415,3 +415,3 @@
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$,
... ...
@@ -420,6 +420,6 @@
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.
Show white space 6 line context
... ...
@@ -17,3 +17,4 @@
17 17

	
18

	
18
nodist_lemon_HEADERS = lemon/config.h	
19
	
19 20
lemon_libemon_la_CXXFLAGS = \
... ...
@@ -59,2 +60,3 @@
59 60
	lemon/bin_heap.h \
61
	lemon/cbc.h \
60 62
	lemon/circulation.h \
... ...
@@ -63,3 +65,2 @@
63 65
	lemon/concept_check.h \
64
	lemon/config.h \
65 66
	lemon/connectivity.h \
Ignore white space 6 line context
... ...
@@ -66,4 +66,4 @@
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

	
... ...
@@ -74,3 +74,3 @@
74 74
    /// concept.
75
    typedef typename Digraph::template ArcMap<Flow> FlowMap;
75
    typedef typename Digraph::template ArcMap<Value> FlowMap;
76 76

	
... ...
@@ -106,3 +106,3 @@
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

	
... ...
@@ -189,4 +189,4 @@
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

	
... ...
@@ -223,3 +223,3 @@
223 223

	
224
    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
224
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
225 225
    ExcessMap* _excess;
... ...
@@ -532,3 +532,3 @@
532 532
        } else {
533
          Flow fc = -(*_excess)[_g.target(e)];
533
          Value fc = -(*_excess)[_g.target(e)];
534 534
          _flow->set(e, fc);
... ...
@@ -565,3 +565,3 @@
565 565
        int mlevel=_node_num;
566
        Flow exc=(*_excess)[act];
566
        Value exc=(*_excess)[act];
567 567

	
... ...
@@ -569,3 +569,3 @@
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;
... ...
@@ -593,3 +593,3 @@
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;
... ...
@@ -663,5 +663,5 @@
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
    ///
... ...
@@ -669,3 +669,3 @@
669 669
    /// using this function.
670
    Flow flow(const Arc& arc) const {
670
    Value flow(const Arc& arc) const {
671 671
      return (*_flow)[arc];
... ...
@@ -752,3 +752,3 @@
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];
... ...
@@ -767,6 +767,6 @@
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)
Ignore white space 6 line context
... ...
@@ -24,3 +24,3 @@
24 24

	
25
#include <lemon/core.h>
25
#include <lemon/config.h>
26 26
#include <lemon/bits/enable_if.h>
Ignore white space 6 line context
... ...
@@ -32,5 +32,2 @@
32 32
#include <lemon/math.h>
33
#include <lemon/maps.h>
34
#include <lemon/circulation.h>
35
#include <lemon/adaptors.h>
36 33

	
... ...
@@ -52,10 +49,15 @@
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
  ///
... ...
@@ -67,3 +69,3 @@
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
... ...
@@ -72,17 +74,6 @@
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

	
... ...
@@ -90,7 +81,76 @@
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
... ...
@@ -133,54 +193,2 @@
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:
... ...
@@ -189,6 +197,2 @@
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;
... ...
@@ -197,3 +201,3 @@
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;
... ...
@@ -215,16 +219,5 @@
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

	
... ...
@@ -232,3 +225,3 @@
232 225
    IntNodeMap _node_id;
233
    ArcVector _arc_ref;
226
    IntArcMap _arc_id;
234 227
    IntVector _source;
... ...
@@ -237,6 +230,8 @@
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;
... ...
@@ -259,3 +254,12 @@
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

	
... ...
@@ -661,21 +665,64 @@
661 665
    NetworkSimplex(const GR& graph) :
662
      _graph(graph),
663
      _plower(NULL), _pupper(NULL), _pcost(NULL),
664
      _psupply(NULL), _pstsup(false), _ptype(GEQ),
665
      _flow_map(NULL), _potential_map(NULL),
666
      _local_flow(false), _local_potential(false),
667
      _node_id(graph)
666
      _graph(graph), _node_id(graph), _arc_id(graph),
667
      INF(std::numeric_limits<Value>::has_infinity ?
668
          std::numeric_limits<Value>::infinity() :
669
          std::numeric_limits<Value>::max())
668 670
    {
669
      LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
670
                   std::numeric_limits<Flow>::is_signed,
671
        "The flow type of NetworkSimplex must be signed integer");
672
      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
673
                   std::numeric_limits<Cost>::is_signed,
674
        "The cost type of NetworkSimplex must be signed integer");
675
    }
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;
676 682

	
677
    /// Destructor.
678
    ~NetworkSimplex() {
679
      if (_local_flow) delete _flow_map;
680
      if (_local_potential) delete _potential_map;
683
      _source.resize(all_arc_num);
684
      _target.resize(all_arc_num);
685

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

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

	
703
      // Copy the graph (store the arcs in a mixed order)
704
      int i = 0;
705
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
706
        _node_id[n] = i;
707
      }
708
      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
709
      i = 0;
710
      for (ArcIt a(_graph); a != INVALID; ++a) {
711
        _arc_id[a] = i;
712
        _source[i] = _node_id[_graph.source(a)];
713
        _target[i] = _node_id[_graph.target(a)];
714
        if ((i += k) >= _arc_num) i = (i % k) + 1;
715
      }
716
      
717
      // Initialize maps
718
      for (int i = 0; i != _node_num; ++i) {
719
        _supply[i] = 0;
720
      }
721
      for (int i = 0; i != _arc_num; ++i) {
722
        _lower[i] = 0;
723
        _upper[i] = INF;
724
        _cost[i] = 1;
725
      }
726
      _have_lower = false;
727
      _stype = GEQ;
681 728
    }
... ...
@@ -691,8 +738,7 @@
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.
... ...
@@ -700,8 +746,7 @@
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
      }
... ...
@@ -713,9 +758,8 @@
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.
... ...
@@ -723,8 +767,6 @@
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
      }
... ...
@@ -733,39 +775,2 @@
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.
... ...
@@ -781,8 +786,6 @@
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
      }
... ...
@@ -799,3 +802,3 @@
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.
... ...
@@ -803,9 +806,6 @@
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
      }
... ...
@@ -822,2 +822,6 @@
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.
... ...
@@ -828,9 +832,8 @@
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;
... ...
@@ -838,13 +841,13 @@
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;
851
    NetworkSimplex& supplyType(SupplyType supply_type) {
852
      _stype = supply_type;
850 853
      return *this;
... ...
@@ -852,37 +855,2 @@
852 855

	
853
    /// \brief Set the flow map.
854
    ///
855
    /// This function sets the flow map.
856
    /// If it is not used before calling \ref run(), an instance 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;
885
      return *this;
886
    }
887
    
888 856
    /// @}
... ...
@@ -898,5 +866,4 @@
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,
... ...
@@ -904,3 +871,3 @@
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();
... ...
@@ -912,2 +879,4 @@
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
    ///
... ...
@@ -916,5 +885,14 @@
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
    }
... ...
@@ -925,5 +903,3 @@
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
    ///
... ...
@@ -932,2 +908,4 @@
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
    ///
... ...
@@ -938,3 +916,3 @@
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();
... ...
@@ -949,3 +927,3 @@
949 927
    ///   ns.reset();
950
    ///   ns.capacityMap(cap).costMap(cost)
928
    ///   ns.upperMap(capacity).costMap(cost)
951 929
    ///     .supplyMap(sup).run();
... ...
@@ -955,19 +933,12 @@
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;
... ...
@@ -987,3 +958,3 @@
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
    ///
... ...
@@ -999,11 +970,8 @@
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
      }
... ...
@@ -1023,14 +991,18 @@
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
    }
... ...
@@ -1044,15 +1016,18 @@
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
    }
... ...
@@ -1065,152 +1040,48 @@
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;
1109
      }
1110
      if (_psupply) {
1111
        int i = 0;
1112
        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1113
          _node_id[n] = i;
1114
          _supply[i] = (*_psupply)[n];
1115
          sum_supply += _supply[i];
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];
1116 1077
        }
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
... ...
@@ -1221,65 +1092,6 @@
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

	
... ...
@@ -1287,2 +1099,4 @@
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;
... ...
@@ -1291,11 +1105,9 @@
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 {
... ...
@@ -1303,3 +1115,3 @@
1303 1115
          _forward[u] = false;
1304
          _pi[u] = art_cost + _pi[_root];
1116
          _pi[u] = ART_COST + _pi[_root];
1305 1117
        }
... ...
@@ -1338,3 +1150,3 @@
1338 1150
      int result = 0;
1339
      Flow d;
1151
      Value d;
1340 1152
      int e;
... ...
@@ -1344,3 +1156,4 @@
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) {
... ...
@@ -1354,3 +1167,4 @@
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) {
... ...
@@ -1376,3 +1190,3 @@
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;
... ...
@@ -1528,3 +1342,3 @@
1528 1342
    // Execute the algorithm
1529
    bool start(PivotRule pivot_rule) {
1343
    ProblemType start(PivotRule pivot_rule) {
1530 1344
      // Select the pivot rule implementation
... ...
@@ -1542,3 +1356,3 @@
1542 1356
      }
1543
      return false;
1357
      return INFEASIBLE; // avoid warning
1544 1358
    }
... ...
@@ -1546,3 +1360,3 @@
1546 1360
    template <typename PivotRuleImpl>
1547
    bool start() {
1361
    ProblemType start() {
1548 1362
      PivotRuleImpl pivot(*this);
... ...
@@ -1553,2 +1367,3 @@
1553 1367
        bool change = findLeavingArc();
1368
        if (delta >= INF) return UNBOUNDED;
1554 1369
        changeFlow(change);
... ...
@@ -1559,20 +1374,33 @@
1559 1374
      }
1560

	
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]);
1566
        }
1567
      } else {
1568
        for (int i = 0; i != _arc_num; ++i) {
1569
          _flow_map->set(_arc_ref[i], _flow[i]);
1375
      
1376
      // Check feasibility
1377
      if (_sum_supply < 0) {
1378
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1379
          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
1570 1380
        }
1571 1381
      }
1572
      // Copy potential values to _potential_map
1573
      for (NodeIt n(_graph); n != INVALID; ++n) {
1574
        _potential_map->set(n, _pi[_node_id[n]]);
1382
      else if (_sum_supply > 0) {
1383
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1384
          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
1385
        }
1386
      }
1387
      else {
1388
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1389
          if (_flow[e] != 0) return INFEASIBLE;
1390
        }
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
    }
Ignore white space 6 line context
... ...
@@ -48,3 +48,3 @@
48 48
    /// \brief The type of the flow values.
49
    typedef typename CapacityMap::Value Flow;
49
    typedef typename CapacityMap::Value Value;
50 50

	
... ...
@@ -54,3 +54,3 @@
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

	
... ...
@@ -86,3 +86,3 @@
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

	
... ...
@@ -127,3 +127,3 @@
127 127
    ///The type of the flow values.
128
    typedef typename Traits::Flow Flow;
128
    typedef typename Traits::Value Value;
129 129

	
... ...
@@ -153,3 +153,3 @@
153 153

	
154
    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
154
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
155 155
    ExcessMap* _excess;
... ...
@@ -472,3 +472,3 @@
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) {
... ...
@@ -521,3 +521,3 @@
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)) {
... ...
@@ -533,3 +533,3 @@
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)) {
... ...
@@ -566,3 +566,3 @@
566 566
        while (num > 0 && n != INVALID) {
567
          Flow excess = (*_excess)[n];
567
          Value excess = (*_excess)[n];
568 568
          int new_level = _level->maxLevel();
... ...
@@ -570,3 +570,3 @@
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;
... ...
@@ -593,3 +593,3 @@
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;
... ...
@@ -639,3 +639,3 @@
639 639
        while (num > 0 && n != INVALID) {
640
          Flow excess = (*_excess)[n];
640
          Value excess = (*_excess)[n];
641 641
          int new_level = _level->maxLevel();
... ...
@@ -643,3 +643,3 @@
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;
... ...
@@ -666,3 +666,3 @@
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;
... ...
@@ -780,3 +780,3 @@
780 780
      while ((n = _level->highestActive()) != INVALID) {
781
        Flow excess = (*_excess)[n];
781
        Value excess = (*_excess)[n];
782 782
        int level = _level->highestActiveLevel();
... ...
@@ -785,3 +785,3 @@
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;
... ...
@@ -808,3 +808,3 @@
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;
... ...
@@ -899,3 +899,3 @@
899 899
    /// using this function.
900
    Flow flowValue() const {
900
    Value flowValue() const {
901 901
      return (*_excess)[_target];
... ...
@@ -903,5 +903,5 @@
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.
... ...
@@ -910,3 +910,3 @@
910 910
    /// using this function.
911
    Flow flow(const Arc& arc) const {
911
    Value flow(const Arc& arc) const {
912 912
      return (*_flow)[arc];
Ignore white space 6 line context
... ...
@@ -23,5 +23,3 @@
23 23

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

	
Ignore white space 6 line context
... ...
@@ -20,2 +20,3 @@
20 20
#include <fstream>
21
#include <limits>
21 22

	
... ...
@@ -35,39 +36,39 @@
35 36
  "@nodes\n"
36
  "label  sup1 sup2 sup3 sup4 sup5\n"
37
  "    1    20   27    0   20   30\n"
38
  "    2    -4    0    0   -8   -3\n"
39
  "    3     0    0    0    0    0\n"
40
  "    4     0    0    0    0    0\n"
41
  "    5     9    0    0    6   11\n"
42
  "    6    -6    0    0   -5   -6\n"
43
  "    7     0    0    0    0    0\n"
44
  "    8     0    0    0    0    3\n"
45
  "    9     3    0    0    0    0\n"
46
  "   10    -2    0    0   -7   -2\n"
47
  "   11     0    0    0  -10    0\n"
48
  "   12   -20  -27    0  -30  -20\n"
49
  "\n"
37
  "label  sup1 sup2 sup3 sup4 sup5 sup6\n"
38
  "    1    20   27    0   30   20   30\n"
39
  "    2    -4    0    0    0   -8   -3\n"
40
  "    3     0    0    0    0    0    0\n"
41
  "    4     0    0    0    0    0    0\n"
42
  "    5     9    0    0    0    6   11\n"
43
  "    6    -6    0    0    0   -5   -6\n"
44
  "    7     0    0    0    0    0    0\n"
45
  "    8     0    0    0    0    0    3\n"
46
  "    9     3    0    0    0    0    0\n"
47
  "   10    -2    0    0    0   -7   -2\n"
48
  "   11     0    0    0    0  -10    0\n"
49
  "   12   -20  -27    0  -30  -30  -20\n"
50
  "\n"                
50 51
  "@arcs\n"
51
  "       cost  cap low1 low2\n"
52
  " 1  2    70   11    0    8\n"
53
  " 1  3   150    3    0    1\n"
54
  " 1  4    80   15    0    2\n"
55
  " 2  8    80   12    0    0\n"
56
  " 3  5   140    5    0    3\n"
57
  " 4  6    60   10    0    1\n"
58
  " 4  7    80    2    0    0\n"
59
  " 4  8   110    3    0    0\n"
60
  " 5  7    60   14    0    0\n"
61
  " 5 11   120   12    0    0\n"
62
  " 6  3     0    3    0    0\n"
63
  " 6  9   140    4    0    0\n"
64
  " 6 10    90    8    0    0\n"
65
  " 7  1    30    5    0    0\n"
66
  " 8 12    60   16    0    4\n"
67
  " 9 12    50    6    0    0\n"
68
  "10 12    70   13    0    5\n"
69
  "10  2   100    7    0    0\n"
70
  "10  7    60   10    0    0\n"
71
  "11 10    20   14    0    6\n"
72
  "12 11    30   10    0    0\n"
52
  "       cost  cap low1 low2 low3\n"
53
  " 1  2    70   11    0    8    8\n"
54
  " 1  3   150    3    0    1    0\n"
55
  " 1  4    80   15    0    2    2\n"
56
  " 2  8    80   12    0    0    0\n"
57
  " 3  5   140    5    0    3    1\n"
58
  " 4  6    60   10    0    1    0\n"
59
  " 4  7    80    2    0    0    0\n"
60
  " 4  8   110    3    0    0    0\n"
61
  " 5  7    60   14    0    0    0\n"
62
  " 5 11   120   12    0    0    0\n"
63
  " 6  3     0    3    0    0    0\n"
64
  " 6  9   140    4    0    0    0\n"
65
  " 6 10    90    8    0    0    0\n"
66
  " 7  1    30    5    0    0   -5\n"
67
  " 8 12    60   16    0    4    3\n"
68
  " 9 12    50    6    0    0    0\n"
69
  "10 12    70   13    0    5    2\n"
70
  "10  2   100    7    0    0    0\n"
71
  "10  7    60   10    0    0   -3\n"
72
  "11 10    20   14    0    6  -20\n"
73
  "12 11    30   10    0    0  -10\n"
73 74
  "\n"
... ...
@@ -78,3 +79,3 @@
78 79

	
79
enum ProblemType {
80
enum SupplyType {
80 81
  EQ,
... ...
@@ -85,3 +86,3 @@
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
... ...
@@ -96,2 +97,3 @@
96 97
      MCF mcf(g);
98
      const MCF& const_mcf = mcf;
97 99

	
... ...
@@ -100,4 +102,2 @@
100 102
             .upperMap(upper)
101
             .capacityMap(upper)
102
             .boundMaps(lower, upper)
103 103
             .costMap(cost)
... ...
@@ -105,19 +105,10 @@
105 105
             .stSupply(n, n, k)
106
             .flowMap(flow)
107
             .potentialMap(pot)
108 106
             .run();
109
      
110
      const MCF& const_mcf = mcf;
111 107

	
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
    }
... ...
@@ -126,9 +117,11 @@
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;
... ...
@@ -137,8 +130,9 @@
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
  };
... ...
@@ -153,3 +147,3 @@
153 147
                const SM& supply, const FM& flow,
154
                ProblemType type = EQ )
148
                SupplyType type = EQ )
155 149
{
... ...
@@ -210,17 +204,21 @@
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);
... ...
@@ -233,7 +231,9 @@
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
  }
... ...
@@ -246,4 +246,4 @@
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());
... ...
@@ -257,2 +257,3 @@
257 257
    .arcMap("low2", l2)
258
    .arcMap("low3", l3)
258 259
    .nodeMap("sup1", s1)
... ...
@@ -262,2 +263,3 @@
262 263
    .nodeMap("sup5", s5)
264
    .nodeMap("sup6", s6)
263 265
    .node("source", v)
... ...
@@ -265,2 +267,42 @@
265 267
    .run();
268
  
269
  // Build a test digraph for testing negative costs
270
  Digraph ngr;
271
  Node n1 = ngr.addNode();
272
  Node n2 = ngr.addNode();
273
  Node n3 = ngr.addNode();
274
  Node n4 = ngr.addNode();
275
  Node n5 = ngr.addNode();
276
  Node n6 = ngr.addNode();
277
  Node n7 = ngr.addNode();
278
  
279
  Arc a1 = ngr.addArc(n1, n2);
280
  Arc a2 = ngr.addArc(n1, n3);
281
  Arc a3 = ngr.addArc(n2, n4);
282
  Arc a4 = ngr.addArc(n3, n4);
283
  Arc a5 = ngr.addArc(n3, n2);
284
  Arc a6 = ngr.addArc(n5, n3);
285
  Arc a7 = ngr.addArc(n5, n6);
286
  Arc a8 = ngr.addArc(n6, n7);
287
  Arc a9 = ngr.addArc(n7, n5);
288
  
289
  Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
290
  ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
291
  Digraph::NodeMap<int> ns(ngr, 0);
292
  
293
  nl2[a7] =  1000;
294
  nl2[a8] = -1000;
295
  
296
  ns[n1] =  100;
297
  ns[n4] = -100;
298
  
299
  nc[a1] =  100;
300
  nc[a2] =   30;
301
  nc[a3] =   20;
302
  nc[a4] =   80;
303
  nc[a5] =   50;
304
  nc[a6] =   10;
305
  nc[a7] =   80;
306
  nc[a8] =   30;
307
  nc[a9] = -120;
266 308

	
... ...
@@ -273,42 +315,56 @@
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
  }
... ...
@@ -318,14 +374,14 @@
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
  }
Ignore white space 6 line context
... ...
@@ -20,5 +20,3 @@
20 20

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

	
Ignore white space 6 line context
... ...
@@ -121,4 +121,4 @@
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';
0 comments (0 inline)