Support negative costs and bounds in NetworkSimplex (#270)
authorPeter Kovacs <kpeter@inf.elte.hu>
Wed, 29 Apr 2009 03:15:24 +0200
changeset 6366c408d864fa1
parent 618 58357e986a08
child 637 756a5ec551c8
Support negative costs and bounds in NetworkSimplex (#270)

* The interface is reworked to support negative costs and bounds.
- ProblemType and problemType() are renamed to
SupplyType and supplyType(), see also #234.
- ProblemType type is introduced similarly to the LP interface.
- 'bool run()' is replaced by 'ProblemType run()' to handle
unbounded problem instances, as well.
- Add INF public member constant similarly to the LP interface.
* Remove capacityMap() and boundMaps(), see also #266.
* Update the problem definition in the MCF module.
* Remove the usage of Circulation (and adaptors) for checking feasibility.
Check feasibility by examining the artifical arcs instead (after solving
the problem).
* Additional check for unbounded negative cycles found during the
algorithm (it is possible now, since negative costs are allowed).
* Fix in the constructor (the value types needn't be integer any more),
see also #254.
* Improve and extend the doc.
* Rework the test file and add test cases for negative costs and bounds.
doc/groups.dox
lemon/network_simplex.h
test/min_cost_flow_test.cc
tools/dimacs-solver.cc
     1.1 --- a/doc/groups.dox	Sun Apr 26 16:36:23 2009 +0100
     1.2 +++ b/doc/groups.dox	Wed Apr 29 03:15:24 2009 +0200
     1.3 @@ -352,17 +352,17 @@
     1.4  minimum total cost from a set of supply nodes to a set of demand nodes
     1.5  in a network with capacity constraints (lower and upper bounds)
     1.6  and arc costs.
     1.7 -Formally, let \f$G=(V,A)\f$ be a digraph,
     1.8 -\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
     1.9 +Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
    1.10 +\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
    1.11  upper bounds for the flow values on the arcs, for which
    1.12 -\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
    1.13 -\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
    1.14 -on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
    1.15 +\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
    1.16 +\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
    1.17 +on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
    1.18  signed supply values of the nodes.
    1.19  If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
    1.20  supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
    1.21  \f$-sup(u)\f$ demand.
    1.22 -A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
    1.23 +A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
    1.24  of the following optimization problem.
    1.25  
    1.26  \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    1.27 @@ -404,7 +404,7 @@
    1.28  
    1.29  The dual solution of the minimum cost flow problem is represented by node 
    1.30  potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
    1.31 -An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
    1.32 +An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
    1.33  is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
    1.34  node potentials the following \e complementary \e slackness optimality
    1.35  conditions hold.
    1.36 @@ -413,15 +413,15 @@
    1.37     - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
    1.38     - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
    1.39     - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
    1.40 - - For all \f$u\in V\f$:
    1.41 + - For all \f$u\in V\f$ nodes:
    1.42     - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
    1.43       then \f$\pi(u)=0\f$.
    1.44   
    1.45  Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
    1.46 -\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
    1.47 +\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
    1.48  \f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
    1.49  
    1.50 -All algorithms provide dual solution (node potentials) as well
    1.51 +All algorithms provide dual solution (node potentials) as well,
    1.52  if an optimal flow is found.
    1.53  
    1.54  LEMON contains several algorithms for solving minimum cost flow problems.
     2.1 --- a/lemon/network_simplex.h	Sun Apr 26 16:36:23 2009 +0100
     2.2 +++ b/lemon/network_simplex.h	Wed Apr 29 03:15:24 2009 +0200
     2.3 @@ -30,9 +30,6 @@
     2.4  
     2.5  #include <lemon/core.h>
     2.6  #include <lemon/math.h>
     2.7 -#include <lemon/maps.h>
     2.8 -#include <lemon/circulation.h>
     2.9 -#include <lemon/adaptors.h>
    2.10  
    2.11  namespace lemon {
    2.12  
    2.13 @@ -50,8 +47,13 @@
    2.14    ///
    2.15    /// In general this class is the fastest implementation available
    2.16    /// in LEMON for the minimum cost flow problem.
    2.17 -  /// Moreover it supports both direction of the supply/demand inequality
    2.18 -  /// constraints. For more information see \ref ProblemType.
    2.19 +  /// Moreover it supports both directions of the supply/demand inequality
    2.20 +  /// constraints. For more information see \ref SupplyType.
    2.21 +  ///
    2.22 +  /// Most of the parameters of the problem (except for the digraph)
    2.23 +  /// can be given using separate functions, and the algorithm can be
    2.24 +  /// executed using the \ref run() function. If some parameters are not
    2.25 +  /// specified, then default values will be used.
    2.26    ///
    2.27    /// \tparam GR The digraph type the algorithm runs on.
    2.28    /// \tparam F The value type used for flow amounts, capacity bounds
    2.29 @@ -88,11 +90,80 @@
    2.30  
    2.31    public:
    2.32  
    2.33 -    /// \brief Enum type for selecting the pivot rule.
    2.34 +    /// \brief Problem type constants for the \c run() function.
    2.35      ///
    2.36 -    /// Enum type for selecting the pivot rule for the \ref run()
    2.37 +    /// Enum type containing the problem type constants that can be
    2.38 +    /// returned by the \ref run() function of the algorithm.
    2.39 +    enum ProblemType {
    2.40 +      /// The problem has no feasible solution (flow).
    2.41 +      INFEASIBLE,
    2.42 +      /// The problem has optimal solution (i.e. it is feasible and
    2.43 +      /// bounded), and the algorithm has found optimal flow and node
    2.44 +      /// potentials (primal and dual solutions).
    2.45 +      OPTIMAL,
    2.46 +      /// The objective function of the problem is unbounded, i.e.
    2.47 +      /// there is a directed cycle having negative total cost and
    2.48 +      /// infinite upper bound.
    2.49 +      UNBOUNDED
    2.50 +    };
    2.51 +    
    2.52 +    /// \brief Constants for selecting the type of the supply constraints.
    2.53 +    ///
    2.54 +    /// Enum type containing constants for selecting the supply type,
    2.55 +    /// i.e. the direction of the inequalities in the supply/demand
    2.56 +    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
    2.57 +    ///
    2.58 +    /// The default supply type is \c GEQ, since this form is supported
    2.59 +    /// by other minimum cost flow algorithms and the \ref Circulation
    2.60 +    /// algorithm, as well.
    2.61 +    /// The \c LEQ problem type can be selected using the \ref supplyType()
    2.62      /// function.
    2.63      ///
    2.64 +    /// Note that the equality form is a special case of both supply types.
    2.65 +    enum SupplyType {
    2.66 +
    2.67 +      /// This option means that there are <em>"greater or equal"</em>
    2.68 +      /// supply/demand constraints in the definition, i.e. the exact
    2.69 +      /// formulation of the problem is the following.
    2.70 +      /**
    2.71 +          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    2.72 +          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
    2.73 +              sup(u) \quad \forall u\in V \f]
    2.74 +          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    2.75 +      */
    2.76 +      /// It means that the total demand must be greater or equal to the 
    2.77 +      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
    2.78 +      /// negative) and all the supplies have to be carried out from 
    2.79 +      /// the supply nodes, but there could be demands that are not 
    2.80 +      /// satisfied.
    2.81 +      GEQ,
    2.82 +      /// It is just an alias for the \c GEQ option.
    2.83 +      CARRY_SUPPLIES = GEQ,
    2.84 +
    2.85 +      /// This option means that there are <em>"less or equal"</em>
    2.86 +      /// supply/demand constraints in the definition, i.e. the exact
    2.87 +      /// formulation of the problem is the following.
    2.88 +      /**
    2.89 +          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    2.90 +          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
    2.91 +              sup(u) \quad \forall u\in V \f]
    2.92 +          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    2.93 +      */
    2.94 +      /// It means that the total demand must be less or equal to the 
    2.95 +      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
    2.96 +      /// positive) and all the demands have to be satisfied, but there
    2.97 +      /// could be supplies that are not carried out from the supply
    2.98 +      /// nodes.
    2.99 +      LEQ,
   2.100 +      /// It is just an alias for the \c LEQ option.
   2.101 +      SATISFY_DEMANDS = LEQ
   2.102 +    };
   2.103 +    
   2.104 +    /// \brief Constants for selecting the pivot rule.
   2.105 +    ///
   2.106 +    /// Enum type containing constants for selecting the pivot rule for
   2.107 +    /// the \ref run() function.
   2.108 +    ///
   2.109      /// \ref NetworkSimplex provides five different pivot rule
   2.110      /// implementations that significantly affect the running time
   2.111      /// of the algorithm.
   2.112 @@ -131,58 +202,6 @@
   2.113        ALTERING_LIST
   2.114      };
   2.115      
   2.116 -    /// \brief Enum type for selecting the problem type.
   2.117 -    ///
   2.118 -    /// Enum type for selecting the problem type, i.e. the direction of
   2.119 -    /// the inequalities in the supply/demand constraints of the
   2.120 -    /// \ref min_cost_flow "minimum cost flow problem".
   2.121 -    ///
   2.122 -    /// The default problem type is \c GEQ, since this form is supported
   2.123 -    /// by other minimum cost flow algorithms and the \ref Circulation
   2.124 -    /// algorithm as well.
   2.125 -    /// The \c LEQ problem type can be selected using the \ref problemType()
   2.126 -    /// function.
   2.127 -    ///
   2.128 -    /// Note that the equality form is a special case of both problem type.
   2.129 -    enum ProblemType {
   2.130 -
   2.131 -      /// This option means that there are "<em>greater or equal</em>"
   2.132 -      /// constraints in the defintion, i.e. the exact formulation of the
   2.133 -      /// problem is the following.
   2.134 -      /**
   2.135 -          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
   2.136 -          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
   2.137 -              sup(u) \quad \forall u\in V \f]
   2.138 -          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
   2.139 -      */
   2.140 -      /// It means that the total demand must be greater or equal to the 
   2.141 -      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
   2.142 -      /// negative) and all the supplies have to be carried out from 
   2.143 -      /// the supply nodes, but there could be demands that are not 
   2.144 -      /// satisfied.
   2.145 -      GEQ,
   2.146 -      /// It is just an alias for the \c GEQ option.
   2.147 -      CARRY_SUPPLIES = GEQ,
   2.148 -
   2.149 -      /// This option means that there are "<em>less or equal</em>"
   2.150 -      /// constraints in the defintion, i.e. the exact formulation of the
   2.151 -      /// problem is the following.
   2.152 -      /**
   2.153 -          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
   2.154 -          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
   2.155 -              sup(u) \quad \forall u\in V \f]
   2.156 -          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
   2.157 -      */
   2.158 -      /// It means that the total demand must be less or equal to the 
   2.159 -      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
   2.160 -      /// positive) and all the demands have to be satisfied, but there
   2.161 -      /// could be supplies that are not carried out from the supply
   2.162 -      /// nodes.
   2.163 -      LEQ,
   2.164 -      /// It is just an alias for the \c LEQ option.
   2.165 -      SATISFY_DEMANDS = LEQ
   2.166 -    };
   2.167 -
   2.168    private:
   2.169  
   2.170      TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   2.171 @@ -220,7 +239,9 @@
   2.172      bool _pstsup;
   2.173      Node _psource, _ptarget;
   2.174      Flow _pstflow;
   2.175 -    ProblemType _ptype;
   2.176 +    SupplyType _stype;
   2.177 +    
   2.178 +    Flow _sum_supply;
   2.179  
   2.180      // Result maps
   2.181      FlowMap *_flow_map;
   2.182 @@ -259,6 +280,15 @@
   2.183      int stem, par_stem, new_stem;
   2.184      Flow delta;
   2.185  
   2.186 +  public:
   2.187 +  
   2.188 +    /// \brief Constant for infinite upper bounds (capacities).
   2.189 +    ///
   2.190 +    /// Constant for infinite upper bounds (capacities).
   2.191 +    /// It is \c std::numeric_limits<Flow>::infinity() if available,
   2.192 +    /// \c std::numeric_limits<Flow>::max() otherwise.
   2.193 +    const Flow INF;
   2.194 +
   2.195    private:
   2.196  
   2.197      // Implementation of the First Eligible pivot rule
   2.198 @@ -661,17 +691,19 @@
   2.199      NetworkSimplex(const GR& graph) :
   2.200        _graph(graph),
   2.201        _plower(NULL), _pupper(NULL), _pcost(NULL),
   2.202 -      _psupply(NULL), _pstsup(false), _ptype(GEQ),
   2.203 +      _psupply(NULL), _pstsup(false), _stype(GEQ),
   2.204        _flow_map(NULL), _potential_map(NULL),
   2.205        _local_flow(false), _local_potential(false),
   2.206 -      _node_id(graph)
   2.207 +      _node_id(graph),
   2.208 +      INF(std::numeric_limits<Flow>::has_infinity ?
   2.209 +          std::numeric_limits<Flow>::infinity() :
   2.210 +          std::numeric_limits<Flow>::max())
   2.211      {
   2.212 -      LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
   2.213 -                   std::numeric_limits<Flow>::is_signed,
   2.214 -        "The flow type of NetworkSimplex must be signed integer");
   2.215 -      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
   2.216 -                   std::numeric_limits<Cost>::is_signed,
   2.217 -        "The cost type of NetworkSimplex must be signed integer");
   2.218 +      // Check the value types
   2.219 +      LEMON_ASSERT(std::numeric_limits<Flow>::is_signed,
   2.220 +        "The flow type of NetworkSimplex must be signed");
   2.221 +      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
   2.222 +        "The cost type of NetworkSimplex must be signed");
   2.223      }
   2.224  
   2.225      /// Destructor.
   2.226 @@ -689,17 +721,16 @@
   2.227      /// \brief Set the lower bounds on the arcs.
   2.228      ///
   2.229      /// This function sets the lower bounds on the arcs.
   2.230 -    /// If neither this function nor \ref boundMaps() is used before
   2.231 -    /// calling \ref run(), the lower bounds will be set to zero
   2.232 -    /// on all arcs.
   2.233 +    /// If it is not used before calling \ref run(), the lower bounds
   2.234 +    /// will be set to zero on all arcs.
   2.235      ///
   2.236      /// \param map An arc map storing the lower bounds.
   2.237      /// Its \c Value type must be convertible to the \c Flow type
   2.238      /// of the algorithm.
   2.239      ///
   2.240      /// \return <tt>(*this)</tt>
   2.241 -    template <typename LOWER>
   2.242 -    NetworkSimplex& lowerMap(const LOWER& map) {
   2.243 +    template <typename LowerMap>
   2.244 +    NetworkSimplex& lowerMap(const LowerMap& map) {
   2.245        delete _plower;
   2.246        _plower = new FlowArcMap(_graph);
   2.247        for (ArcIt a(_graph); a != INVALID; ++a) {
   2.248 @@ -711,18 +742,17 @@
   2.249      /// \brief Set the upper bounds (capacities) on the arcs.
   2.250      ///
   2.251      /// This function sets the upper bounds (capacities) on the arcs.
   2.252 -    /// If none of the functions \ref upperMap(), \ref capacityMap()
   2.253 -    /// and \ref boundMaps() is used before calling \ref run(),
   2.254 -    /// the upper bounds (capacities) will be set to
   2.255 -    /// \c std::numeric_limits<Flow>::max() on all arcs.
   2.256 +    /// If it is not used before calling \ref run(), the upper bounds
   2.257 +    /// will be set to \ref INF on all arcs (i.e. the flow value will be
   2.258 +    /// unbounded from above on each arc).
   2.259      ///
   2.260      /// \param map An arc map storing the upper bounds.
   2.261      /// Its \c Value type must be convertible to the \c Flow type
   2.262      /// of the algorithm.
   2.263      ///
   2.264      /// \return <tt>(*this)</tt>
   2.265 -    template<typename UPPER>
   2.266 -    NetworkSimplex& upperMap(const UPPER& map) {
   2.267 +    template<typename UpperMap>
   2.268 +    NetworkSimplex& upperMap(const UpperMap& map) {
   2.269        delete _pupper;
   2.270        _pupper = new FlowArcMap(_graph);
   2.271        for (ArcIt a(_graph); a != INVALID; ++a) {
   2.272 @@ -731,43 +761,6 @@
   2.273        return *this;
   2.274      }
   2.275  
   2.276 -    /// \brief Set the upper bounds (capacities) on the arcs.
   2.277 -    ///
   2.278 -    /// This function sets the upper bounds (capacities) on the arcs.
   2.279 -    /// It is just an alias for \ref upperMap().
   2.280 -    ///
   2.281 -    /// \return <tt>(*this)</tt>
   2.282 -    template<typename CAP>
   2.283 -    NetworkSimplex& capacityMap(const CAP& map) {
   2.284 -      return upperMap(map);
   2.285 -    }
   2.286 -
   2.287 -    /// \brief Set the lower and upper bounds on the arcs.
   2.288 -    ///
   2.289 -    /// This function sets the lower and upper bounds on the arcs.
   2.290 -    /// If neither this function nor \ref lowerMap() is used before
   2.291 -    /// calling \ref run(), the lower bounds will be set to zero
   2.292 -    /// on all arcs.
   2.293 -    /// If none of the functions \ref upperMap(), \ref capacityMap()
   2.294 -    /// and \ref boundMaps() is used before calling \ref run(),
   2.295 -    /// the upper bounds (capacities) will be set to
   2.296 -    /// \c std::numeric_limits<Flow>::max() on all arcs.
   2.297 -    ///
   2.298 -    /// \param lower An arc map storing the lower bounds.
   2.299 -    /// \param upper An arc map storing the upper bounds.
   2.300 -    ///
   2.301 -    /// The \c Value type of the maps must be convertible to the
   2.302 -    /// \c Flow type of the algorithm.
   2.303 -    ///
   2.304 -    /// \note This function is just a shortcut of calling \ref lowerMap()
   2.305 -    /// and \ref upperMap() separately.
   2.306 -    ///
   2.307 -    /// \return <tt>(*this)</tt>
   2.308 -    template <typename LOWER, typename UPPER>
   2.309 -    NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
   2.310 -      return lowerMap(lower).upperMap(upper);
   2.311 -    }
   2.312 -
   2.313      /// \brief Set the costs of the arcs.
   2.314      ///
   2.315      /// This function sets the costs of the arcs.
   2.316 @@ -779,8 +772,8 @@
   2.317      /// of the algorithm.
   2.318      ///
   2.319      /// \return <tt>(*this)</tt>
   2.320 -    template<typename COST>
   2.321 -    NetworkSimplex& costMap(const COST& map) {
   2.322 +    template<typename CostMap>
   2.323 +    NetworkSimplex& costMap(const CostMap& map) {
   2.324        delete _pcost;
   2.325        _pcost = new CostArcMap(_graph);
   2.326        for (ArcIt a(_graph); a != INVALID; ++a) {
   2.327 @@ -801,8 +794,8 @@
   2.328      /// of the algorithm.
   2.329      ///
   2.330      /// \return <tt>(*this)</tt>
   2.331 -    template<typename SUP>
   2.332 -    NetworkSimplex& supplyMap(const SUP& map) {
   2.333 +    template<typename SupplyMap>
   2.334 +    NetworkSimplex& supplyMap(const SupplyMap& map) {
   2.335        delete _psupply;
   2.336        _pstsup = false;
   2.337        _psupply = new FlowNodeMap(_graph);
   2.338 @@ -820,6 +813,10 @@
   2.339      /// calling \ref run(), the supply of each node will be set to zero.
   2.340      /// (It makes sense only if non-zero lower bounds are given.)
   2.341      ///
   2.342 +    /// Using this function has the same effect as using \ref supplyMap()
   2.343 +    /// with such a map in which \c k is assigned to \c s, \c -k is
   2.344 +    /// assigned to \c t and all other nodes have zero supply value.
   2.345 +    ///
   2.346      /// \param s The source node.
   2.347      /// \param t The target node.
   2.348      /// \param k The required amount of flow from node \c s to node \c t
   2.349 @@ -836,17 +833,17 @@
   2.350        return *this;
   2.351      }
   2.352      
   2.353 -    /// \brief Set the problem type.
   2.354 +    /// \brief Set the type of the supply constraints.
   2.355      ///
   2.356 -    /// This function sets the problem type for the algorithm.
   2.357 -    /// If it is not used before calling \ref run(), the \ref GEQ problem
   2.358 +    /// This function sets the type of the supply/demand constraints.
   2.359 +    /// If it is not used before calling \ref run(), the \ref GEQ supply
   2.360      /// type will be used.
   2.361      ///
   2.362 -    /// For more information see \ref ProblemType.
   2.363 +    /// For more information see \ref SupplyType.
   2.364      ///
   2.365      /// \return <tt>(*this)</tt>
   2.366 -    NetworkSimplex& problemType(ProblemType problem_type) {
   2.367 -      _ptype = problem_type;
   2.368 +    NetworkSimplex& supplyType(SupplyType supply_type) {
   2.369 +      _stype = supply_type;
   2.370        return *this;
   2.371      }
   2.372  
   2.373 @@ -896,13 +893,12 @@
   2.374      ///
   2.375      /// This function runs the algorithm.
   2.376      /// The paramters can be specified using functions \ref lowerMap(),
   2.377 -    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
   2.378 -    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
   2.379 -    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
   2.380 +    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 
   2.381 +    /// \ref supplyType(), \ref flowMap() and \ref potentialMap().
   2.382      /// For example,
   2.383      /// \code
   2.384      ///   NetworkSimplex<ListDigraph> ns(graph);
   2.385 -    ///   ns.boundMaps(lower, upper).costMap(cost)
   2.386 +    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
   2.387      ///     .supplyMap(sup).run();
   2.388      /// \endcode
   2.389      ///
   2.390 @@ -914,17 +910,25 @@
   2.391      /// \param pivot_rule The pivot rule that will be used during the
   2.392      /// algorithm. For more information see \ref PivotRule.
   2.393      ///
   2.394 -    /// \return \c true if a feasible flow can be found.
   2.395 -    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
   2.396 -      return init() && start(pivot_rule);
   2.397 +    /// \return \c INFEASIBLE if no feasible flow exists,
   2.398 +    /// \n \c OPTIMAL if the problem has optimal solution
   2.399 +    /// (i.e. it is feasible and bounded), and the algorithm has found
   2.400 +    /// optimal flow and node potentials (primal and dual solutions),
   2.401 +    /// \n \c UNBOUNDED if the objective function of the problem is
   2.402 +    /// unbounded, i.e. there is a directed cycle having negative total
   2.403 +    /// cost and infinite upper bound.
   2.404 +    ///
   2.405 +    /// \see ProblemType, PivotRule
   2.406 +    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
   2.407 +      if (!init()) return INFEASIBLE;
   2.408 +      return start(pivot_rule);
   2.409      }
   2.410  
   2.411      /// \brief Reset all the parameters that have been given before.
   2.412      ///
   2.413      /// This function resets all the paramaters that have been given
   2.414      /// before using functions \ref lowerMap(), \ref upperMap(),
   2.415 -    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
   2.416 -    /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
   2.417 +    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(),
   2.418      /// \ref flowMap() and \ref potentialMap().
   2.419      ///
   2.420      /// It is useful for multiple run() calls. If this function is not
   2.421 @@ -936,7 +940,7 @@
   2.422      ///   NetworkSimplex<ListDigraph> ns(graph);
   2.423      ///
   2.424      ///   // First run
   2.425 -    ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
   2.426 +    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
   2.427      ///     .supplyMap(sup).run();
   2.428      ///
   2.429      ///   // Run again with modified cost map (reset() is not called,
   2.430 @@ -947,7 +951,7 @@
   2.431      ///   // Run again from scratch using reset()
   2.432      ///   // (the lower bounds will be set to zero on all arcs)
   2.433      ///   ns.reset();
   2.434 -    ///   ns.capacityMap(cap).costMap(cost)
   2.435 +    ///   ns.upperMap(capacity).costMap(cost)
   2.436      ///     .supplyMap(sup).run();
   2.437      /// \endcode
   2.438      ///
   2.439 @@ -962,7 +966,7 @@
   2.440        _pcost = NULL;
   2.441        _psupply = NULL;
   2.442        _pstsup = false;
   2.443 -      _ptype = GEQ;
   2.444 +      _stype = GEQ;
   2.445        if (_local_flow) delete _flow_map;
   2.446        if (_local_potential) delete _potential_map;
   2.447        _flow_map = NULL;
   2.448 @@ -985,7 +989,7 @@
   2.449      /// \brief Return the total cost of the found flow.
   2.450      ///
   2.451      /// This function returns the total cost of the found flow.
   2.452 -    /// The complexity of the function is O(e).
   2.453 +    /// Its complexity is O(e).
   2.454      ///
   2.455      /// \note The return type of the function can be specified as a
   2.456      /// template parameter. For example,
   2.457 @@ -997,9 +1001,9 @@
   2.458      /// function.
   2.459      ///
   2.460      /// \pre \ref run() must be called before using this function.
   2.461 -    template <typename Num>
   2.462 -    Num totalCost() const {
   2.463 -      Num c = 0;
   2.464 +    template <typename Value>
   2.465 +    Value totalCost() const {
   2.466 +      Value c = 0;
   2.467        if (_pcost) {
   2.468          for (ArcIt e(_graph); e != INVALID; ++e)
   2.469            c += (*_flow_map)[e] * (*_pcost)[e];
   2.470 @@ -1050,7 +1054,7 @@
   2.471      ///
   2.472      /// This function returns a const reference to a node map storing
   2.473      /// the found potentials, which form the dual solution of the
   2.474 -    /// \ref min_cost_flow "minimum cost flow" problem.
   2.475 +    /// \ref min_cost_flow "minimum cost flow problem".
   2.476      ///
   2.477      /// \pre \ref run() must be called before using this function.
   2.478      const PotentialMap& potentialMap() const {
   2.479 @@ -1101,7 +1105,7 @@
   2.480  
   2.481        // Initialize node related data
   2.482        bool valid_supply = true;
   2.483 -      Flow sum_supply = 0;
   2.484 +      _sum_supply = 0;
   2.485        if (!_pstsup && !_psupply) {
   2.486          _pstsup = true;
   2.487          _psource = _ptarget = NodeIt(_graph);
   2.488 @@ -1112,10 +1116,10 @@
   2.489          for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   2.490            _node_id[n] = i;
   2.491            _supply[i] = (*_psupply)[n];
   2.492 -          sum_supply += _supply[i];
   2.493 +          _sum_supply += _supply[i];
   2.494          }
   2.495 -        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
   2.496 -                       (_ptype == LEQ && sum_supply >= 0);
   2.497 +        valid_supply = (_stype == GEQ && _sum_supply <= 0) ||
   2.498 +                       (_stype == LEQ && _sum_supply >= 0);
   2.499        } else {
   2.500          int i = 0;
   2.501          for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   2.502 @@ -1127,92 +1131,18 @@
   2.503        }
   2.504        if (!valid_supply) return false;
   2.505  
   2.506 -      // Infinite capacity value
   2.507 -      Flow inf_cap =
   2.508 -        std::numeric_limits<Flow>::has_infinity ?
   2.509 -        std::numeric_limits<Flow>::infinity() :
   2.510 -        std::numeric_limits<Flow>::max();
   2.511 -
   2.512        // Initialize artifical cost
   2.513 -      Cost art_cost;
   2.514 +      Cost ART_COST;
   2.515        if (std::numeric_limits<Cost>::is_exact) {
   2.516 -        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
   2.517 +        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
   2.518        } else {
   2.519 -        art_cost = std::numeric_limits<Cost>::min();
   2.520 +        ART_COST = std::numeric_limits<Cost>::min();
   2.521          for (int i = 0; i != _arc_num; ++i) {
   2.522 -          if (_cost[i] > art_cost) art_cost = _cost[i];
   2.523 +          if (_cost[i] > ART_COST) ART_COST = _cost[i];
   2.524          }
   2.525 -        art_cost = (art_cost + 1) * _node_num;
   2.526 +        ART_COST = (ART_COST + 1) * _node_num;
   2.527        }
   2.528  
   2.529 -      // Run Circulation to check if a feasible solution exists
   2.530 -      typedef ConstMap<Arc, Flow> ConstArcMap;
   2.531 -      ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
   2.532 -      FlowNodeMap *csup = NULL;
   2.533 -      bool local_csup = false;
   2.534 -      if (_psupply) {
   2.535 -        csup = _psupply;
   2.536 -      } else {
   2.537 -        csup = new FlowNodeMap(_graph, 0);
   2.538 -        (*csup)[_psource] =  _pstflow;
   2.539 -        (*csup)[_ptarget] = -_pstflow;
   2.540 -        local_csup = true;
   2.541 -      }
   2.542 -      bool circ_result = false;
   2.543 -      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
   2.544 -        // GEQ problem type
   2.545 -        if (_plower) {
   2.546 -          if (_pupper) {
   2.547 -            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
   2.548 -              circ(_graph, *_plower, *_pupper, *csup);
   2.549 -            circ_result = circ.run();
   2.550 -          } else {
   2.551 -            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
   2.552 -              circ(_graph, *_plower, inf_arc_map, *csup);
   2.553 -            circ_result = circ.run();
   2.554 -          }
   2.555 -        } else {
   2.556 -          if (_pupper) {
   2.557 -            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
   2.558 -              circ(_graph, zero_arc_map, *_pupper, *csup);
   2.559 -            circ_result = circ.run();
   2.560 -          } else {
   2.561 -            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
   2.562 -              circ(_graph, zero_arc_map, inf_arc_map, *csup);
   2.563 -            circ_result = circ.run();
   2.564 -          }
   2.565 -        }
   2.566 -      } else {
   2.567 -        // LEQ problem type
   2.568 -        typedef ReverseDigraph<const GR> RevGraph;
   2.569 -        typedef NegMap<FlowNodeMap> NegNodeMap;
   2.570 -        RevGraph rgraph(_graph);
   2.571 -        NegNodeMap neg_csup(*csup);
   2.572 -        if (_plower) {
   2.573 -          if (_pupper) {
   2.574 -            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
   2.575 -              circ(rgraph, *_plower, *_pupper, neg_csup);
   2.576 -            circ_result = circ.run();
   2.577 -          } else {
   2.578 -            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
   2.579 -              circ(rgraph, *_plower, inf_arc_map, neg_csup);
   2.580 -            circ_result = circ.run();
   2.581 -          }
   2.582 -        } else {
   2.583 -          if (_pupper) {
   2.584 -            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
   2.585 -              circ(rgraph, zero_arc_map, *_pupper, neg_csup);
   2.586 -            circ_result = circ.run();
   2.587 -          } else {
   2.588 -            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
   2.589 -              circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
   2.590 -            circ_result = circ.run();
   2.591 -          }
   2.592 -        }
   2.593 -      }
   2.594 -      if (local_csup) delete csup;
   2.595 -      if (!circ_result) return false;
   2.596 -
   2.597        // Set data for the artificial root node
   2.598        _root = _node_num;
   2.599        _parent[_root] = -1;
   2.600 @@ -1221,11 +1151,11 @@
   2.601        _rev_thread[0] = _root;
   2.602        _succ_num[_root] = all_node_num;
   2.603        _last_succ[_root] = _root - 1;
   2.604 -      _supply[_root] = -sum_supply;
   2.605 -      if (sum_supply < 0) {
   2.606 -        _pi[_root] = -art_cost;
   2.607 +      _supply[_root] = -_sum_supply;
   2.608 +      if (_sum_supply < 0) {
   2.609 +        _pi[_root] = -ART_COST;
   2.610        } else {
   2.611 -        _pi[_root] = art_cost;
   2.612 +        _pi[_root] = ART_COST;
   2.613        }
   2.614  
   2.615        // Store the arcs in a mixed order
   2.616 @@ -1260,7 +1190,7 @@
   2.617              _cap[i] = (*_pupper)[_arc_ref[i]];
   2.618          } else {
   2.619            for (int i = 0; i != _arc_num; ++i)
   2.620 -            _cap[i] = inf_cap;
   2.621 +            _cap[i] = INF;
   2.622          }
   2.623          if (_pcost) {
   2.624            for (int i = 0; i != _arc_num; ++i)
   2.625 @@ -1275,8 +1205,17 @@
   2.626        if (_plower) {
   2.627          for (int i = 0; i != _arc_num; ++i) {
   2.628            Flow c = (*_plower)[_arc_ref[i]];
   2.629 -          if (c != 0) {
   2.630 -            _cap[i] -= c;
   2.631 +          if (c > 0) {
   2.632 +            if (_cap[i] < INF) _cap[i] -= c;
   2.633 +            _supply[_source[i]] -= c;
   2.634 +            _supply[_target[i]] += c;
   2.635 +          }
   2.636 +          else if (c < 0) {
   2.637 +            if (_cap[i] < INF + c) {
   2.638 +              _cap[i] -= c;
   2.639 +            } else {
   2.640 +              _cap[i] = INF;
   2.641 +            }
   2.642              _supply[_source[i]] -= c;
   2.643              _supply[_target[i]] += c;
   2.644            }
   2.645 @@ -1291,17 +1230,17 @@
   2.646          _last_succ[u] = u;
   2.647          _parent[u] = _root;
   2.648          _pred[u] = e;
   2.649 -        _cost[e] = art_cost;
   2.650 -        _cap[e] = inf_cap;
   2.651 +        _cost[e] = ART_COST;
   2.652 +        _cap[e] = INF;
   2.653          _state[e] = STATE_TREE;
   2.654 -        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
   2.655 +        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
   2.656            _flow[e] = _supply[u];
   2.657            _forward[u] = true;
   2.658 -          _pi[u] = -art_cost + _pi[_root];
   2.659 +          _pi[u] = -ART_COST + _pi[_root];
   2.660          } else {
   2.661            _flow[e] = -_supply[u];
   2.662            _forward[u] = false;
   2.663 -          _pi[u] = art_cost + _pi[_root];
   2.664 +          _pi[u] = ART_COST + _pi[_root];
   2.665          }
   2.666        }
   2.667  
   2.668 @@ -1342,7 +1281,8 @@
   2.669        // Search the cycle along the path form the first node to the root
   2.670        for (int u = first; u != join; u = _parent[u]) {
   2.671          e = _pred[u];
   2.672 -        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
   2.673 +        d = _forward[u] ?
   2.674 +          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
   2.675          if (d < delta) {
   2.676            delta = d;
   2.677            u_out = u;
   2.678 @@ -1352,7 +1292,8 @@
   2.679        // Search the cycle along the path form the second node to the root
   2.680        for (int u = second; u != join; u = _parent[u]) {
   2.681          e = _pred[u];
   2.682 -        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
   2.683 +        d = _forward[u] ? 
   2.684 +          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
   2.685          if (d <= delta) {
   2.686            delta = d;
   2.687            u_out = u;
   2.688 @@ -1526,7 +1467,7 @@
   2.689      }
   2.690  
   2.691      // Execute the algorithm
   2.692 -    bool start(PivotRule pivot_rule) {
   2.693 +    ProblemType start(PivotRule pivot_rule) {
   2.694        // Select the pivot rule implementation
   2.695        switch (pivot_rule) {
   2.696          case FIRST_ELIGIBLE:
   2.697 @@ -1540,23 +1481,41 @@
   2.698          case ALTERING_LIST:
   2.699            return start<AlteringListPivotRule>();
   2.700        }
   2.701 -      return false;
   2.702 +      return INFEASIBLE; // avoid warning
   2.703      }
   2.704  
   2.705      template <typename PivotRuleImpl>
   2.706 -    bool start() {
   2.707 +    ProblemType start() {
   2.708        PivotRuleImpl pivot(*this);
   2.709  
   2.710        // Execute the Network Simplex algorithm
   2.711        while (pivot.findEnteringArc()) {
   2.712          findJoinNode();
   2.713          bool change = findLeavingArc();
   2.714 +        if (delta >= INF) return UNBOUNDED;
   2.715          changeFlow(change);
   2.716          if (change) {
   2.717            updateTreeStructure();
   2.718            updatePotential();
   2.719          }
   2.720        }
   2.721 +      
   2.722 +      // Check feasibility
   2.723 +      if (_sum_supply < 0) {
   2.724 +        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
   2.725 +          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
   2.726 +        }
   2.727 +      }
   2.728 +      else if (_sum_supply > 0) {
   2.729 +        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
   2.730 +          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
   2.731 +        }
   2.732 +      }
   2.733 +      else {
   2.734 +        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
   2.735 +          if (_flow[e] != 0) return INFEASIBLE;
   2.736 +        }
   2.737 +      }
   2.738  
   2.739        // Copy flow values to _flow_map
   2.740        if (_plower) {
   2.741 @@ -1574,7 +1533,7 @@
   2.742          _potential_map->set(n, _pi[_node_id[n]]);
   2.743        }
   2.744  
   2.745 -      return true;
   2.746 +      return OPTIMAL;
   2.747      }
   2.748  
   2.749    }; //class NetworkSimplex
     3.1 --- a/test/min_cost_flow_test.cc	Sun Apr 26 16:36:23 2009 +0100
     3.2 +++ b/test/min_cost_flow_test.cc	Wed Apr 29 03:15:24 2009 +0200
     3.3 @@ -18,6 +18,7 @@
     3.4  
     3.5  #include <iostream>
     3.6  #include <fstream>
     3.7 +#include <limits>
     3.8  
     3.9  #include <lemon/list_graph.h>
    3.10  #include <lemon/lgf_reader.h>
    3.11 @@ -33,50 +34,50 @@
    3.12  
    3.13  char test_lgf[] =
    3.14    "@nodes\n"
    3.15 -  "label  sup1 sup2 sup3 sup4 sup5\n"
    3.16 -  "    1    20   27    0   20   30\n"
    3.17 -  "    2    -4    0    0   -8   -3\n"
    3.18 -  "    3     0    0    0    0    0\n"
    3.19 -  "    4     0    0    0    0    0\n"
    3.20 -  "    5     9    0    0    6   11\n"
    3.21 -  "    6    -6    0    0   -5   -6\n"
    3.22 -  "    7     0    0    0    0    0\n"
    3.23 -  "    8     0    0    0    0    3\n"
    3.24 -  "    9     3    0    0    0    0\n"
    3.25 -  "   10    -2    0    0   -7   -2\n"
    3.26 -  "   11     0    0    0  -10    0\n"
    3.27 -  "   12   -20  -27    0  -30  -20\n"
    3.28 -  "\n"
    3.29 +  "label  sup1 sup2 sup3 sup4 sup5 sup6\n"
    3.30 +  "    1    20   27    0   30   20   30\n"
    3.31 +  "    2    -4    0    0    0   -8   -3\n"
    3.32 +  "    3     0    0    0    0    0    0\n"
    3.33 +  "    4     0    0    0    0    0    0\n"
    3.34 +  "    5     9    0    0    0    6   11\n"
    3.35 +  "    6    -6    0    0    0   -5   -6\n"
    3.36 +  "    7     0    0    0    0    0    0\n"
    3.37 +  "    8     0    0    0    0    0    3\n"
    3.38 +  "    9     3    0    0    0    0    0\n"
    3.39 +  "   10    -2    0    0    0   -7   -2\n"
    3.40 +  "   11     0    0    0    0  -10    0\n"
    3.41 +  "   12   -20  -27    0  -30  -30  -20\n"
    3.42 +  "\n"                
    3.43    "@arcs\n"
    3.44 -  "       cost  cap low1 low2\n"
    3.45 -  " 1  2    70   11    0    8\n"
    3.46 -  " 1  3   150    3    0    1\n"
    3.47 -  " 1  4    80   15    0    2\n"
    3.48 -  " 2  8    80   12    0    0\n"
    3.49 -  " 3  5   140    5    0    3\n"
    3.50 -  " 4  6    60   10    0    1\n"
    3.51 -  " 4  7    80    2    0    0\n"
    3.52 -  " 4  8   110    3    0    0\n"
    3.53 -  " 5  7    60   14    0    0\n"
    3.54 -  " 5 11   120   12    0    0\n"
    3.55 -  " 6  3     0    3    0    0\n"
    3.56 -  " 6  9   140    4    0    0\n"
    3.57 -  " 6 10    90    8    0    0\n"
    3.58 -  " 7  1    30    5    0    0\n"
    3.59 -  " 8 12    60   16    0    4\n"
    3.60 -  " 9 12    50    6    0    0\n"
    3.61 -  "10 12    70   13    0    5\n"
    3.62 -  "10  2   100    7    0    0\n"
    3.63 -  "10  7    60   10    0    0\n"
    3.64 -  "11 10    20   14    0    6\n"
    3.65 -  "12 11    30   10    0    0\n"
    3.66 +  "       cost  cap low1 low2 low3\n"
    3.67 +  " 1  2    70   11    0    8    8\n"
    3.68 +  " 1  3   150    3    0    1    0\n"
    3.69 +  " 1  4    80   15    0    2    2\n"
    3.70 +  " 2  8    80   12    0    0    0\n"
    3.71 +  " 3  5   140    5    0    3    1\n"
    3.72 +  " 4  6    60   10    0    1    0\n"
    3.73 +  " 4  7    80    2    0    0    0\n"
    3.74 +  " 4  8   110    3    0    0    0\n"
    3.75 +  " 5  7    60   14    0    0    0\n"
    3.76 +  " 5 11   120   12    0    0    0\n"
    3.77 +  " 6  3     0    3    0    0    0\n"
    3.78 +  " 6  9   140    4    0    0    0\n"
    3.79 +  " 6 10    90    8    0    0    0\n"
    3.80 +  " 7  1    30    5    0    0   -5\n"
    3.81 +  " 8 12    60   16    0    4    3\n"
    3.82 +  " 9 12    50    6    0    0    0\n"
    3.83 +  "10 12    70   13    0    5    2\n"
    3.84 +  "10  2   100    7    0    0    0\n"
    3.85 +  "10  7    60   10    0    0   -3\n"
    3.86 +  "11 10    20   14    0    6  -20\n"
    3.87 +  "12 11    30   10    0    0  -10\n"
    3.88    "\n"
    3.89    "@attributes\n"
    3.90    "source 1\n"
    3.91    "target 12\n";
    3.92  
    3.93  
    3.94 -enum ProblemType {
    3.95 +enum SupplyType {
    3.96    EQ,
    3.97    GEQ,
    3.98    LEQ
    3.99 @@ -98,8 +99,6 @@
   3.100        b = mcf.reset()
   3.101               .lowerMap(lower)
   3.102               .upperMap(upper)
   3.103 -             .capacityMap(upper)
   3.104 -             .boundMaps(lower, upper)
   3.105               .costMap(cost)
   3.106               .supplyMap(sup)
   3.107               .stSupply(n, n, k)
   3.108 @@ -112,10 +111,12 @@
   3.109        const typename MCF::FlowMap &fm = const_mcf.flowMap();
   3.110        const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
   3.111  
   3.112 -      v = const_mcf.totalCost();
   3.113 +      c = const_mcf.totalCost();
   3.114        double x = const_mcf.template totalCost<double>();
   3.115        v = const_mcf.flow(a);
   3.116 -      v = const_mcf.potential(n);
   3.117 +      c = const_mcf.potential(n);
   3.118 +      
   3.119 +      v = const_mcf.INF;
   3.120  
   3.121        ignore_unused_variable_warning(fm);
   3.122        ignore_unused_variable_warning(pm);
   3.123 @@ -137,6 +138,7 @@
   3.124      const Arc &a;
   3.125      const Flow &k;
   3.126      Flow v;
   3.127 +    Cost c;
   3.128      bool b;
   3.129  
   3.130      typename MCF::FlowMap &flow;
   3.131 @@ -151,7 +153,7 @@
   3.132             typename SM, typename FM >
   3.133  bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
   3.134                  const SM& supply, const FM& flow,
   3.135 -                ProblemType type = EQ )
   3.136 +                SupplyType type = EQ )
   3.137  {
   3.138    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   3.139  
   3.140 @@ -208,16 +210,17 @@
   3.141  // Run a minimum cost flow algorithm and check the results
   3.142  template < typename MCF, typename GR,
   3.143             typename LM, typename UM,
   3.144 -           typename CM, typename SM >
   3.145 -void checkMcf( const MCF& mcf, bool mcf_result,
   3.146 +           typename CM, typename SM,
   3.147 +           typename PT >
   3.148 +void checkMcf( const MCF& mcf, PT mcf_result,
   3.149                 const GR& gr, const LM& lower, const UM& upper,
   3.150                 const CM& cost, const SM& supply,
   3.151 -               bool result, typename CM::Value total,
   3.152 +               PT result, bool optimal, typename CM::Value total,
   3.153                 const std::string &test_id = "",
   3.154 -               ProblemType type = EQ )
   3.155 +               SupplyType type = EQ )
   3.156  {
   3.157    check(mcf_result == result, "Wrong result " + test_id);
   3.158 -  if (result) {
   3.159 +  if (optimal) {
   3.160      check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
   3.161            "The flow is not feasible " + test_id);
   3.162      check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
   3.163 @@ -244,8 +247,8 @@
   3.164  
   3.165    // Read the test digraph
   3.166    Digraph gr;
   3.167 -  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
   3.168 -  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
   3.169 +  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
   3.170 +  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
   3.171    ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
   3.172    Node v, w;
   3.173  
   3.174 @@ -255,14 +258,56 @@
   3.175      .arcMap("cap", u)
   3.176      .arcMap("low1", l1)
   3.177      .arcMap("low2", l2)
   3.178 +    .arcMap("low3", l3)
   3.179      .nodeMap("sup1", s1)
   3.180      .nodeMap("sup2", s2)
   3.181      .nodeMap("sup3", s3)
   3.182      .nodeMap("sup4", s4)
   3.183      .nodeMap("sup5", s5)
   3.184 +    .nodeMap("sup6", s6)
   3.185      .node("source", v)
   3.186      .node("target", w)
   3.187      .run();
   3.188 +  
   3.189 +  // Build a test digraph for testing negative costs
   3.190 +  Digraph ngr;
   3.191 +  Node n1 = ngr.addNode();
   3.192 +  Node n2 = ngr.addNode();
   3.193 +  Node n3 = ngr.addNode();
   3.194 +  Node n4 = ngr.addNode();
   3.195 +  Node n5 = ngr.addNode();
   3.196 +  Node n6 = ngr.addNode();
   3.197 +  Node n7 = ngr.addNode();
   3.198 +  
   3.199 +  Arc a1 = ngr.addArc(n1, n2);
   3.200 +  Arc a2 = ngr.addArc(n1, n3);
   3.201 +  Arc a3 = ngr.addArc(n2, n4);
   3.202 +  Arc a4 = ngr.addArc(n3, n4);
   3.203 +  Arc a5 = ngr.addArc(n3, n2);
   3.204 +  Arc a6 = ngr.addArc(n5, n3);
   3.205 +  Arc a7 = ngr.addArc(n5, n6);
   3.206 +  Arc a8 = ngr.addArc(n6, n7);
   3.207 +  Arc a9 = ngr.addArc(n7, n5);
   3.208 +  
   3.209 +  Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
   3.210 +  ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
   3.211 +  Digraph::NodeMap<int> ns(ngr, 0);
   3.212 +  
   3.213 +  nl2[a7] =  1000;
   3.214 +  nl2[a8] = -1000;
   3.215 +  
   3.216 +  ns[n1] =  100;
   3.217 +  ns[n4] = -100;
   3.218 +  
   3.219 +  nc[a1] =  100;
   3.220 +  nc[a2] =   30;
   3.221 +  nc[a3] =   20;
   3.222 +  nc[a4] =   80;
   3.223 +  nc[a5] =   50;
   3.224 +  nc[a6] =   10;
   3.225 +  nc[a7] =   80;
   3.226 +  nc[a8] =   30;
   3.227 +  nc[a9] = -120;
   3.228  
   3.229    // A. Test NetworkSimplex with the default pivot rule
   3.230    {
   3.231 @@ -271,63 +316,77 @@
   3.232      // Check the equality form
   3.233      mcf.upperMap(u).costMap(c);
   3.234      checkMcf(mcf, mcf.supplyMap(s1).run(),
   3.235 -             gr, l1, u, c, s1, true,  5240, "#A1");
   3.236 +             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
   3.237      checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   3.238 -             gr, l1, u, c, s2, true,  7620, "#A2");
   3.239 +             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
   3.240      mcf.lowerMap(l2);
   3.241      checkMcf(mcf, mcf.supplyMap(s1).run(),
   3.242 -             gr, l2, u, c, s1, true,  5970, "#A3");
   3.243 +             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
   3.244      checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
   3.245 -             gr, l2, u, c, s2, true,  8010, "#A4");
   3.246 +             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
   3.247      mcf.reset();
   3.248      checkMcf(mcf, mcf.supplyMap(s1).run(),
   3.249 -             gr, l1, cu, cc, s1, true,  74, "#A5");
   3.250 +             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
   3.251      checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
   3.252 -             gr, l2, cu, cc, s2, true,  94, "#A6");
   3.253 +             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
   3.254      mcf.reset();
   3.255      checkMcf(mcf, mcf.run(),
   3.256 -             gr, l1, cu, cc, s3, true,   0, "#A7");
   3.257 -    checkMcf(mcf, mcf.boundMaps(l2, u).run(),
   3.258 -             gr, l2, u, cc, s3, false,   0, "#A8");
   3.259 +             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
   3.260 +    checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
   3.261 +             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
   3.262 +    mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
   3.263 +    checkMcf(mcf, mcf.run(),
   3.264 +             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
   3.265  
   3.266      // Check the GEQ form
   3.267 -    mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
   3.268 +    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
   3.269      checkMcf(mcf, mcf.run(),
   3.270 -             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
   3.271 -    mcf.problemType(mcf.GEQ);
   3.272 +             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
   3.273 +    mcf.supplyType(mcf.GEQ);
   3.274      checkMcf(mcf, mcf.lowerMap(l2).run(),
   3.275 -             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
   3.276 -    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
   3.277 +             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
   3.278 +    mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
   3.279      checkMcf(mcf, mcf.run(),
   3.280 -             gr, l2, u, c, s5, false,    0, "#A11", GEQ);
   3.281 +             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
   3.282  
   3.283      // Check the LEQ form
   3.284 -    mcf.reset().problemType(mcf.LEQ);
   3.285 -    mcf.upperMap(u).costMap(c).supplyMap(s5);
   3.286 +    mcf.reset().supplyType(mcf.LEQ);
   3.287 +    mcf.upperMap(u).costMap(c).supplyMap(s6);
   3.288      checkMcf(mcf, mcf.run(),
   3.289 -             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
   3.290 +             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
   3.291      checkMcf(mcf, mcf.lowerMap(l2).run(),
   3.292 -             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
   3.293 -    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
   3.294 +             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
   3.295 +    mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
   3.296      checkMcf(mcf, mcf.run(),
   3.297 -             gr, l2, u, c, s4, false,    0, "#A14", LEQ);
   3.298 +             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
   3.299 +
   3.300 +    // Check negative costs
   3.301 +    NetworkSimplex<Digraph> nmcf(ngr);
   3.302 +    nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
   3.303 +    checkMcf(nmcf, nmcf.run(),
   3.304 +      ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
   3.305 +    checkMcf(nmcf, nmcf.upperMap(nu2).run(),
   3.306 +      ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
   3.307 +    nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
   3.308 +    checkMcf(nmcf, nmcf.run(),
   3.309 +      ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
   3.310    }
   3.311  
   3.312    // B. Test NetworkSimplex with each pivot rule
   3.313    {
   3.314      NetworkSimplex<Digraph> mcf(gr);
   3.315 -    mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
   3.316 +    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
   3.317  
   3.318      checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
   3.319 -             gr, l2, u, c, s1, true,  5970, "#B1");
   3.320 +             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
   3.321      checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
   3.322 -             gr, l2, u, c, s1, true,  5970, "#B2");
   3.323 +             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
   3.324      checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
   3.325 -             gr, l2, u, c, s1, true,  5970, "#B3");
   3.326 +             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
   3.327      checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
   3.328 -             gr, l2, u, c, s1, true,  5970, "#B4");
   3.329 +             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
   3.330      checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
   3.331 -             gr, l2, u, c, s1, true,  5970, "#B5");
   3.332 +             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
   3.333    }
   3.334  
   3.335    return 0;
     4.1 --- a/tools/dimacs-solver.cc	Sun Apr 26 16:36:23 2009 +0100
     4.2 +++ b/tools/dimacs-solver.cc	Wed Apr 29 03:15:24 2009 +0200
     4.3 @@ -119,8 +119,8 @@
     4.4  
     4.5    ti.restart();
     4.6    NetworkSimplex<Digraph, Value> ns(g);
     4.7 -  ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
     4.8 -  if (sum_sup > 0) ns.problemType(ns.LEQ);
     4.9 +  ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup);
    4.10 +  if (sum_sup > 0) ns.supplyType(ns.LEQ);
    4.11    if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';
    4.12    ti.restart();
    4.13    bool res = ns.run();