Support >= and <= constraints in NetworkSimplex (#219, #234)
authorPeter Kovacs <kpeter@inf.elte.hu>
Fri, 17 Apr 2009 18:04:36 +0200
changeset 601e6927fe719e6
parent 600 6ac5d9ae1d3d
child 602 dacc2cee2b4c
Support >= and <= constraints in NetworkSimplex (#219, #234)

By default the same inequality constraints are supported as by
Circulation (the GEQ form), but the LEQ form can also be selected
using the problemType() function.

The documentation of the min. cost flow module is reworked and
extended with important notes and explanations about the different
variants of the problem and about the dual solution and optimality
conditions.
doc/groups.dox
lemon/network_simplex.h
test/min_cost_flow_test.cc
     1.1 --- a/doc/groups.dox	Fri Apr 03 18:59:15 2009 +0200
     1.2 +++ b/doc/groups.dox	Fri Apr 17 18:04:36 2009 +0200
     1.3 @@ -317,15 +317,15 @@
     1.4  
     1.5  The \e maximum \e flow \e problem is to find a flow of maximum value between
     1.6  a single source and a single target. Formally, there is a \f$G=(V,A)\f$
     1.7 -digraph, a \f$cap:A\rightarrow\mathbf{R}^+_0\f$ capacity function and
     1.8 +digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and
     1.9  \f$s, t \in V\f$ source and target nodes.
    1.10 -A maximum flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of the
    1.11 +A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the
    1.12  following optimization problem.
    1.13  
    1.14 -\f[ \max\sum_{a\in\delta_{out}(s)}f(a) - \sum_{a\in\delta_{in}(s)}f(a) \f]
    1.15 -\f[ \sum_{a\in\delta_{out}(v)} f(a) = \sum_{a\in\delta_{in}(v)} f(a)
    1.16 -    \qquad \forall v\in V\setminus\{s,t\} \f]
    1.17 -\f[ 0 \leq f(a) \leq cap(a) \qquad \forall a\in A \f]
    1.18 +\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f]
    1.19 +\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu)
    1.20 +    \quad \forall u\in V\setminus\{s,t\} \f]
    1.21 +\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f]
    1.22  
    1.23  LEMON contains several algorithms for solving maximum flow problems:
    1.24  - \ref EdmondsKarp Edmonds-Karp algorithm.
    1.25 @@ -345,35 +345,103 @@
    1.26  
    1.27  \brief Algorithms for finding minimum cost flows and circulations.
    1.28  
    1.29 -This group describes the algorithms for finding minimum cost flows and
    1.30 +This group contains the algorithms for finding minimum cost flows and
    1.31  circulations.
    1.32  
    1.33  The \e minimum \e cost \e flow \e problem is to find a feasible flow of
    1.34  minimum total cost from a set of supply nodes to a set of demand nodes
    1.35 -in a network with capacity constraints and arc costs.
    1.36 +in a network with capacity constraints (lower and upper bounds)
    1.37 +and arc costs.
    1.38  Formally, let \f$G=(V,A)\f$ be a digraph,
    1.39  \f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
    1.40 -upper bounds for the flow values on the arcs,
    1.41 +upper bounds for the flow values on the arcs, for which
    1.42 +\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
    1.43  \f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
    1.44 -on the arcs, and
    1.45 -\f$supply: V\rightarrow\mathbf{Z}\f$ denotes the supply/demand values
    1.46 -of the nodes.
    1.47 -A minimum cost flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of
    1.48 -the following optimization problem.
    1.49 +on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
    1.50 +signed supply values of the nodes.
    1.51 +If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
    1.52 +supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
    1.53 +\f$-sup(u)\f$ demand.
    1.54 +A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
    1.55 +of the following optimization problem.
    1.56  
    1.57 -\f[ \min\sum_{a\in A} f(a) cost(a) \f]
    1.58 -\f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a) =
    1.59 -    supply(v) \qquad \forall v\in V \f]
    1.60 -\f[ lower(a) \leq f(a) \leq upper(a) \qquad \forall a\in A \f]
    1.61 +\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    1.62 +\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
    1.63 +    sup(u) \quad \forall u\in V \f]
    1.64 +\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    1.65  
    1.66 -LEMON contains several algorithms for solving minimum cost flow problems:
    1.67 - - \ref CycleCanceling Cycle-canceling algorithms.
    1.68 - - \ref CapacityScaling Successive shortest path algorithm with optional
    1.69 +The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
    1.70 +zero or negative in order to have a feasible solution (since the sum
    1.71 +of the expressions on the left-hand side of the inequalities is zero).
    1.72 +It means that the total demand must be greater or equal to the total
    1.73 +supply and all the supplies have to be carried out from the supply nodes,
    1.74 +but there could be demands that are not satisfied.
    1.75 +If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
    1.76 +constraints have to be satisfied with equality, i.e. all demands
    1.77 +have to be satisfied and all supplies have to be used.
    1.78 +
    1.79 +If you need the opposite inequalities in the supply/demand constraints
    1.80 +(i.e. the total demand is less than the total supply and all the demands
    1.81 +have to be satisfied while there could be supplies that are not used),
    1.82 +then you could easily transform the problem to the above form by reversing
    1.83 +the direction of the arcs and taking the negative of the supply values
    1.84 +(e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
    1.85 +However \ref NetworkSimplex algorithm also supports this form directly
    1.86 +for the sake of convenience.
    1.87 +
    1.88 +A feasible solution for this problem can be found using \ref Circulation.
    1.89 +
    1.90 +Note that the above formulation is actually more general than the usual
    1.91 +definition of the minimum cost flow problem, in which strict equalities
    1.92 +are required in the supply/demand contraints, i.e.
    1.93 +
    1.94 +\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
    1.95 +    sup(u) \quad \forall u\in V. \f]
    1.96 +
    1.97 +However if the sum of the supply values is zero, then these two problems
    1.98 +are equivalent. So if you need the equality form, you have to ensure this
    1.99 +additional contraint for the algorithms.
   1.100 +
   1.101 +The dual solution of the minimum cost flow problem is represented by node 
   1.102 +potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
   1.103 +An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
   1.104 +is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
   1.105 +node potentials the following \e complementary \e slackness optimality
   1.106 +conditions hold.
   1.107 +
   1.108 + - For all \f$uv\in A\f$ arcs:
   1.109 +   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
   1.110 +   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
   1.111 +   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
   1.112 + - For all \f$u\in V\f$:
   1.113 +   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
   1.114 +     then \f$\pi(u)=0\f$.
   1.115 + 
   1.116 +Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
   1.117 +\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
   1.118 +\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
   1.119 +
   1.120 +All algorithms provide dual solution (node potentials) as well
   1.121 +if an optimal flow is found.
   1.122 +
   1.123 +LEMON contains several algorithms for solving minimum cost flow problems.
   1.124 + - \ref NetworkSimplex Primal Network Simplex algorithm with various
   1.125 +   pivot strategies.
   1.126 + - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
   1.127 +   cost scaling.
   1.128 + - \ref CapacityScaling Successive Shortest %Path algorithm with optional
   1.129     capacity scaling.
   1.130 - - \ref CostScaling Push-relabel and augment-relabel algorithms based on
   1.131 -   cost scaling.
   1.132 - - \ref NetworkSimplex Primal network simplex algorithm with various
   1.133 -   pivot strategies.
   1.134 + - \ref CancelAndTighten The Cancel and Tighten algorithm.
   1.135 + - \ref CycleCanceling Cycle-Canceling algorithms.
   1.136 +
   1.137 +Most of these implementations support the general inequality form of the
   1.138 +minimum cost flow problem, but CancelAndTighten and CycleCanceling
   1.139 +only support the equality form due to the primal method they use.
   1.140 +
   1.141 +In general NetworkSimplex is the most efficient implementation,
   1.142 +but in special cases other algorithms could be faster.
   1.143 +For example, if the total supply and/or capacities are rather small,
   1.144 +CapacityScaling is usually the fastest algorithm (without effective scaling).
   1.145  */
   1.146  
   1.147  /**
     2.1 --- a/lemon/network_simplex.h	Fri Apr 03 18:59:15 2009 +0200
     2.2 +++ b/lemon/network_simplex.h	Fri Apr 17 18:04:36 2009 +0200
     2.3 @@ -30,6 +30,9 @@
     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 @@ -47,6 +50,8 @@
    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    ///
    2.20    /// \tparam GR The digraph type the algorithm runs on.
    2.21    /// \tparam F The value type used for flow amounts, capacity bounds
    2.22 @@ -58,7 +63,8 @@
    2.23    /// be integer.
    2.24    ///
    2.25    /// \note %NetworkSimplex provides five different pivot rule
    2.26 -  /// implementations. For more information see \ref PivotRule.
    2.27 +  /// implementations, from which the most efficient one is used
    2.28 +  /// by default. For more information see \ref PivotRule.
    2.29    template <typename GR, typename F = int, typename C = F>
    2.30    class NetworkSimplex
    2.31    {
    2.32 @@ -68,10 +74,17 @@
    2.33      typedef F Flow;
    2.34      /// The cost type of the algorithm
    2.35      typedef C Cost;
    2.36 +#ifdef DOXYGEN
    2.37 +    /// The type of the flow map
    2.38 +    typedef GR::ArcMap<Flow> FlowMap;
    2.39 +    /// The type of the potential map
    2.40 +    typedef GR::NodeMap<Cost> PotentialMap;
    2.41 +#else
    2.42      /// The type of the flow map
    2.43      typedef typename GR::template ArcMap<Flow> FlowMap;
    2.44      /// The type of the potential map
    2.45      typedef typename GR::template NodeMap<Cost> PotentialMap;
    2.46 +#endif
    2.47  
    2.48    public:
    2.49  
    2.50 @@ -117,6 +130,58 @@
    2.51        /// candidate list and extends this list in every iteration.
    2.52        ALTERING_LIST
    2.53      };
    2.54 +    
    2.55 +    /// \brief Enum type for selecting the problem type.
    2.56 +    ///
    2.57 +    /// Enum type for selecting the problem type, i.e. the direction of
    2.58 +    /// the inequalities in the supply/demand constraints of the
    2.59 +    /// \ref min_cost_flow "minimum cost flow problem".
    2.60 +    ///
    2.61 +    /// The default problem type is \c GEQ, since this form is supported
    2.62 +    /// by other minimum cost flow algorithms and the \ref Circulation
    2.63 +    /// algorithm as well.
    2.64 +    /// The \c LEQ problem type can be selected using the \ref problemType()
    2.65 +    /// function.
    2.66 +    ///
    2.67 +    /// Note that the equality form is a special case of both problem type.
    2.68 +    enum ProblemType {
    2.69 +
    2.70 +      /// This option means that there are "<em>greater or equal</em>"
    2.71 +      /// constraints in the defintion, i.e. the exact formulation of the
    2.72 +      /// problem is the following.
    2.73 +      /**
    2.74 +          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    2.75 +          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
    2.76 +              sup(u) \quad \forall u\in V \f]
    2.77 +          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    2.78 +      */
    2.79 +      /// It means that the total demand must be greater or equal to the 
    2.80 +      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
    2.81 +      /// negative) and all the supplies have to be carried out from 
    2.82 +      /// the supply nodes, but there could be demands that are not 
    2.83 +      /// satisfied.
    2.84 +      GEQ,
    2.85 +      /// It is just an alias for the \c GEQ option.
    2.86 +      CARRY_SUPPLIES = GEQ,
    2.87 +
    2.88 +      /// This option means that there are "<em>less or equal</em>"
    2.89 +      /// constraints in the defintion, i.e. the exact formulation of the
    2.90 +      /// problem is the following.
    2.91 +      /**
    2.92 +          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    2.93 +          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
    2.94 +              sup(u) \quad \forall u\in V \f]
    2.95 +          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    2.96 +      */
    2.97 +      /// It means that the total demand must be less or equal to the 
    2.98 +      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
    2.99 +      /// positive) and all the demands have to be satisfied, but there
   2.100 +      /// could be supplies that are not carried out from the supply
   2.101 +      /// nodes.
   2.102 +      LEQ,
   2.103 +      /// It is just an alias for the \c LEQ option.
   2.104 +      SATISFY_DEMANDS = LEQ
   2.105 +    };
   2.106  
   2.107    private:
   2.108  
   2.109 @@ -155,6 +220,7 @@
   2.110      bool _pstsup;
   2.111      Node _psource, _ptarget;
   2.112      Flow _pstflow;
   2.113 +    ProblemType _ptype;
   2.114  
   2.115      // Result maps
   2.116      FlowMap *_flow_map;
   2.117 @@ -586,13 +652,13 @@
   2.118  
   2.119      /// \brief Constructor.
   2.120      ///
   2.121 -    /// Constructor.
   2.122 +    /// The constructor of the class.
   2.123      ///
   2.124      /// \param graph The digraph the algorithm runs on.
   2.125      NetworkSimplex(const GR& graph) :
   2.126        _graph(graph),
   2.127        _plower(NULL), _pupper(NULL), _pcost(NULL),
   2.128 -      _psupply(NULL), _pstsup(false),
   2.129 +      _psupply(NULL), _pstsup(false), _ptype(GEQ),
   2.130        _flow_map(NULL), _potential_map(NULL),
   2.131        _local_flow(false), _local_potential(false),
   2.132        _node_id(graph)
   2.133 @@ -611,6 +677,12 @@
   2.134        if (_local_potential) delete _potential_map;
   2.135      }
   2.136  
   2.137 +    /// \name Parameters
   2.138 +    /// The parameters of the algorithm can be specified using these
   2.139 +    /// functions.
   2.140 +
   2.141 +    /// @{
   2.142 +
   2.143      /// \brief Set the lower bounds on the arcs.
   2.144      ///
   2.145      /// This function sets the lower bounds on the arcs.
   2.146 @@ -760,6 +832,20 @@
   2.147        _pstflow = k;
   2.148        return *this;
   2.149      }
   2.150 +    
   2.151 +    /// \brief Set the problem type.
   2.152 +    ///
   2.153 +    /// This function sets the problem type for the algorithm.
   2.154 +    /// If it is not used before calling \ref run(), the \ref GEQ problem
   2.155 +    /// type will be used.
   2.156 +    ///
   2.157 +    /// For more information see \ref ProblemType.
   2.158 +    ///
   2.159 +    /// \return <tt>(*this)</tt>
   2.160 +    NetworkSimplex& problemType(ProblemType problem_type) {
   2.161 +      _ptype = problem_type;
   2.162 +      return *this;
   2.163 +    }
   2.164  
   2.165      /// \brief Set the flow map.
   2.166      ///
   2.167 @@ -795,6 +881,8 @@
   2.168        _potential_map = &map;
   2.169        return *this;
   2.170      }
   2.171 +    
   2.172 +    /// @}
   2.173  
   2.174      /// \name Execution Control
   2.175      /// The algorithm can be executed using \ref run().
   2.176 @@ -804,10 +892,11 @@
   2.177      /// \brief Run the algorithm.
   2.178      ///
   2.179      /// This function runs the algorithm.
   2.180 -    /// The paramters can be specified using \ref lowerMap(),
   2.181 +    /// The paramters can be specified using functions \ref lowerMap(),
   2.182      /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
   2.183 -    /// \ref costMap(), \ref supplyMap() and \ref stSupply()
   2.184 -    /// functions. For example,
   2.185 +    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
   2.186 +    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
   2.187 +    /// For example,
   2.188      /// \code
   2.189      ///   NetworkSimplex<ListDigraph> ns(graph);
   2.190      ///   ns.boundMaps(lower, upper).costMap(cost)
   2.191 @@ -830,9 +919,10 @@
   2.192      /// \brief Reset all the parameters that have been given before.
   2.193      ///
   2.194      /// This function resets all the paramaters that have been given
   2.195 -    /// using \ref lowerMap(), \ref upperMap(), \ref capacityMap(),
   2.196 -    /// \ref boundMaps(), \ref costMap(), \ref supplyMap() and
   2.197 -    /// \ref stSupply() functions before.
   2.198 +    /// before using functions \ref lowerMap(), \ref upperMap(),
   2.199 +    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
   2.200 +    /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
   2.201 +    /// \ref flowMap() and \ref potentialMap().
   2.202      ///
   2.203      /// It is useful for multiple run() calls. If this function is not
   2.204      /// used, all the parameters given before are kept for the next
   2.205 @@ -869,6 +959,14 @@
   2.206        _pcost = NULL;
   2.207        _psupply = NULL;
   2.208        _pstsup = false;
   2.209 +      _ptype = GEQ;
   2.210 +      if (_local_flow) delete _flow_map;
   2.211 +      if (_local_potential) delete _potential_map;
   2.212 +      _flow_map = NULL;
   2.213 +      _potential_map = NULL;
   2.214 +      _local_flow = false;
   2.215 +      _local_potential = false;
   2.216 +
   2.217        return *this;
   2.218      }
   2.219  
   2.220 @@ -1000,20 +1098,21 @@
   2.221  
   2.222        // Initialize node related data
   2.223        bool valid_supply = true;
   2.224 +      Flow sum_supply = 0;
   2.225        if (!_pstsup && !_psupply) {
   2.226          _pstsup = true;
   2.227          _psource = _ptarget = NodeIt(_graph);
   2.228          _pstflow = 0;
   2.229        }
   2.230        if (_psupply) {
   2.231 -        Flow sum = 0;
   2.232          int i = 0;
   2.233          for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   2.234            _node_id[n] = i;
   2.235            _supply[i] = (*_psupply)[n];
   2.236 -          sum += _supply[i];
   2.237 +          sum_supply += _supply[i];
   2.238          }
   2.239 -        valid_supply = (sum == 0);
   2.240 +        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
   2.241 +                       (_ptype == LEQ && sum_supply >= 0);
   2.242        } else {
   2.243          int i = 0;
   2.244          for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
   2.245 @@ -1021,10 +1120,95 @@
   2.246            _supply[i] = 0;
   2.247          }
   2.248          _supply[_node_id[_psource]] =  _pstflow;
   2.249 -        _supply[_node_id[_ptarget]]   = -_pstflow;
   2.250 +        _supply[_node_id[_ptarget]] = -_pstflow;
   2.251        }
   2.252        if (!valid_supply) return false;
   2.253  
   2.254 +      // Infinite capacity value
   2.255 +      Flow inf_cap =
   2.256 +        std::numeric_limits<Flow>::has_infinity ?
   2.257 +        std::numeric_limits<Flow>::infinity() :
   2.258 +        std::numeric_limits<Flow>::max();
   2.259 +
   2.260 +      // Initialize artifical cost
   2.261 +      Cost art_cost;
   2.262 +      if (std::numeric_limits<Cost>::is_exact) {
   2.263 +        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
   2.264 +      } else {
   2.265 +        art_cost = std::numeric_limits<Cost>::min();
   2.266 +        for (int i = 0; i != _arc_num; ++i) {
   2.267 +          if (_cost[i] > art_cost) art_cost = _cost[i];
   2.268 +        }
   2.269 +        art_cost = (art_cost + 1) * _node_num;
   2.270 +      }
   2.271 +
   2.272 +      // Run Circulation to check if a feasible solution exists
   2.273 +      typedef ConstMap<Arc, Flow> ConstArcMap;
   2.274 +      FlowNodeMap *csup = NULL;
   2.275 +      bool local_csup = false;
   2.276 +      if (_psupply) {
   2.277 +        csup = _psupply;
   2.278 +      } else {
   2.279 +        csup = new FlowNodeMap(_graph, 0);
   2.280 +        (*csup)[_psource] =  _pstflow;
   2.281 +        (*csup)[_ptarget] = -_pstflow;
   2.282 +        local_csup = true;
   2.283 +      }
   2.284 +      bool circ_result = false;
   2.285 +      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
   2.286 +        // GEQ problem type
   2.287 +        if (_plower) {
   2.288 +          if (_pupper) {
   2.289 +            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
   2.290 +              circ(_graph, *_plower, *_pupper, *csup);
   2.291 +            circ_result = circ.run();
   2.292 +          } else {
   2.293 +            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
   2.294 +              circ(_graph, *_plower, ConstArcMap(inf_cap), *csup);
   2.295 +            circ_result = circ.run();
   2.296 +          }
   2.297 +        } else {
   2.298 +          if (_pupper) {
   2.299 +            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
   2.300 +              circ(_graph, ConstArcMap(0), *_pupper, *csup);
   2.301 +            circ_result = circ.run();
   2.302 +          } else {
   2.303 +            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
   2.304 +              circ(_graph, ConstArcMap(0), ConstArcMap(inf_cap), *csup);
   2.305 +            circ_result = circ.run();
   2.306 +          }
   2.307 +        }
   2.308 +      } else {
   2.309 +        // LEQ problem type
   2.310 +        typedef ReverseDigraph<const GR> RevGraph;
   2.311 +        typedef NegMap<FlowNodeMap> NegNodeMap;
   2.312 +        RevGraph rgraph(_graph);
   2.313 +        NegNodeMap neg_csup(*csup);
   2.314 +        if (_plower) {
   2.315 +          if (_pupper) {
   2.316 +            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
   2.317 +              circ(rgraph, *_plower, *_pupper, neg_csup);
   2.318 +            circ_result = circ.run();
   2.319 +          } else {
   2.320 +            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
   2.321 +              circ(rgraph, *_plower, ConstArcMap(inf_cap), neg_csup);
   2.322 +            circ_result = circ.run();
   2.323 +          }
   2.324 +        } else {
   2.325 +          if (_pupper) {
   2.326 +            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
   2.327 +              circ(rgraph, ConstArcMap(0), *_pupper, neg_csup);
   2.328 +            circ_result = circ.run();
   2.329 +          } else {
   2.330 +            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
   2.331 +              circ(rgraph, ConstArcMap(0), ConstArcMap(inf_cap), neg_csup);
   2.332 +            circ_result = circ.run();
   2.333 +          }
   2.334 +        }
   2.335 +      }
   2.336 +      if (local_csup) delete csup;
   2.337 +      if (!circ_result) return false;
   2.338 +
   2.339        // Set data for the artificial root node
   2.340        _root = _node_num;
   2.341        _parent[_root] = -1;
   2.342 @@ -1033,8 +1217,12 @@
   2.343        _rev_thread[0] = _root;
   2.344        _succ_num[_root] = all_node_num;
   2.345        _last_succ[_root] = _root - 1;
   2.346 -      _supply[_root] = 0;
   2.347 -      _pi[_root] = 0;
   2.348 +      _supply[_root] = -sum_supply;
   2.349 +      if (sum_supply < 0) {
   2.350 +        _pi[_root] = -art_cost;
   2.351 +      } else {
   2.352 +        _pi[_root] = art_cost;
   2.353 +      }
   2.354  
   2.355        // Store the arcs in a mixed order
   2.356        int k = std::max(int(sqrt(_arc_num)), 10);
   2.357 @@ -1045,10 +1233,6 @@
   2.358        }
   2.359  
   2.360        // Initialize arc maps
   2.361 -      Flow inf_cap =
   2.362 -        std::numeric_limits<Flow>::has_infinity ?
   2.363 -        std::numeric_limits<Flow>::infinity() :
   2.364 -        std::numeric_limits<Flow>::max();
   2.365        if (_pupper && _pcost) {
   2.366          for (int i = 0; i != _arc_num; ++i) {
   2.367            Arc e = _arc_ref[i];
   2.368 @@ -1083,18 +1267,6 @@
   2.369          }
   2.370        }
   2.371        
   2.372 -      // Initialize artifical cost
   2.373 -      Cost art_cost;
   2.374 -      if (std::numeric_limits<Cost>::is_exact) {
   2.375 -        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
   2.376 -      } else {
   2.377 -        art_cost = std::numeric_limits<Cost>::min();
   2.378 -        for (int i = 0; i != _arc_num; ++i) {
   2.379 -          if (_cost[i] > art_cost) art_cost = _cost[i];
   2.380 -        }
   2.381 -        art_cost = (art_cost + 1) * _node_num;
   2.382 -      }
   2.383 -
   2.384        // Remove non-zero lower bounds
   2.385        if (_plower) {
   2.386          for (int i = 0; i != _arc_num; ++i) {
   2.387 @@ -1118,14 +1290,14 @@
   2.388          _cost[e] = art_cost;
   2.389          _cap[e] = inf_cap;
   2.390          _state[e] = STATE_TREE;
   2.391 -        if (_supply[u] >= 0) {
   2.392 +        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
   2.393            _flow[e] = _supply[u];
   2.394            _forward[u] = true;
   2.395 -          _pi[u] = -art_cost;
   2.396 +          _pi[u] = -art_cost + _pi[_root];
   2.397          } else {
   2.398            _flow[e] = -_supply[u];
   2.399            _forward[u] = false;
   2.400 -          _pi[u] = art_cost;
   2.401 +          _pi[u] = art_cost + _pi[_root];
   2.402          }
   2.403        }
   2.404  
   2.405 @@ -1382,11 +1554,6 @@
   2.406          }
   2.407        }
   2.408  
   2.409 -      // Check if the flow amount equals zero on all the artificial arcs
   2.410 -      for (int e = _arc_num; e != _arc_num + _node_num; ++e) {
   2.411 -        if (_flow[e] > 0) return false;
   2.412 -      }
   2.413 -
   2.414        // Copy flow values to _flow_map
   2.415        if (_plower) {
   2.416          for (int i = 0; i != _arc_num; ++i) {
     3.1 --- a/test/min_cost_flow_test.cc	Fri Apr 03 18:59:15 2009 +0200
     3.2 +++ b/test/min_cost_flow_test.cc	Fri Apr 17 18:04:36 2009 +0200
     3.3 @@ -33,19 +33,19 @@
     3.4  
     3.5  char test_lgf[] =
     3.6    "@nodes\n"
     3.7 -  "label  sup1 sup2 sup3\n"
     3.8 -  "    1    20   27    0\n"
     3.9 -  "    2    -4    0    0\n"
    3.10 -  "    3     0    0    0\n"
    3.11 -  "    4     0    0    0\n"
    3.12 -  "    5     9    0    0\n"
    3.13 -  "    6    -6    0    0\n"
    3.14 -  "    7     0    0    0\n"
    3.15 -  "    8     0    0    0\n"
    3.16 -  "    9     3    0    0\n"
    3.17 -  "   10    -2    0    0\n"
    3.18 -  "   11     0    0    0\n"
    3.19 -  "   12   -20  -27    0\n"
    3.20 +  "label  sup1 sup2 sup3 sup4 sup5\n"
    3.21 +  "    1    20   27    0   20   30\n"
    3.22 +  "    2    -4    0    0   -8   -3\n"
    3.23 +  "    3     0    0    0    0    0\n"
    3.24 +  "    4     0    0    0    0    0\n"
    3.25 +  "    5     9    0    0    6   11\n"
    3.26 +  "    6    -6    0    0   -5   -6\n"
    3.27 +  "    7     0    0    0    0    0\n"
    3.28 +  "    8     0    0    0    0    3\n"
    3.29 +  "    9     3    0    0    0    0\n"
    3.30 +  "   10    -2    0    0   -7   -2\n"
    3.31 +  "   11     0    0    0  -10    0\n"
    3.32 +  "   12   -20  -27    0  -30  -20\n"
    3.33    "\n"
    3.34    "@arcs\n"
    3.35    "       cost  cap low1 low2\n"
    3.36 @@ -76,6 +76,12 @@
    3.37    "target 12\n";
    3.38  
    3.39  
    3.40 +enum ProblemType {
    3.41 +  EQ,
    3.42 +  GEQ,
    3.43 +  LEQ
    3.44 +};
    3.45 +
    3.46  // Check the interface of an MCF algorithm
    3.47  template <typename GR, typename Flow, typename Cost>
    3.48  class McfClassConcept
    3.49 @@ -97,17 +103,19 @@
    3.50               .costMap(cost)
    3.51               .supplyMap(sup)
    3.52               .stSupply(n, n, k)
    3.53 +             .flowMap(flow)
    3.54 +             .potentialMap(pot)
    3.55               .run();
    3.56 +      
    3.57 +      const MCF& const_mcf = mcf;
    3.58  
    3.59 -      const typename MCF::FlowMap &fm = mcf.flowMap();
    3.60 -      const typename MCF::PotentialMap &pm = mcf.potentialMap();
    3.61 +      const typename MCF::FlowMap &fm = const_mcf.flowMap();
    3.62 +      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
    3.63  
    3.64 -      v = mcf.totalCost();
    3.65 -      double x = mcf.template totalCost<double>();
    3.66 -      v = mcf.flow(a);
    3.67 -      v = mcf.potential(n);
    3.68 -      mcf.flowMap(flow);
    3.69 -      mcf.potentialMap(pot);
    3.70 +      v = const_mcf.totalCost();
    3.71 +      double x = const_mcf.template totalCost<double>();
    3.72 +      v = const_mcf.flow(a);
    3.73 +      v = const_mcf.potential(n);
    3.74  
    3.75        ignore_unused_variable_warning(fm);
    3.76        ignore_unused_variable_warning(pm);
    3.77 @@ -142,7 +150,8 @@
    3.78  template < typename GR, typename LM, typename UM,
    3.79             typename SM, typename FM >
    3.80  bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
    3.81 -                const SM& supply, const FM& flow )
    3.82 +                const SM& supply, const FM& flow,
    3.83 +                ProblemType type = EQ )
    3.84  {
    3.85    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    3.86  
    3.87 @@ -156,7 +165,10 @@
    3.88        sum += flow[e];
    3.89      for (InArcIt e(gr, n); e != INVALID; ++e)
    3.90        sum -= flow[e];
    3.91 -    if (sum != supply[n]) return false;
    3.92 +    bool b = (type ==  EQ && sum == supply[n]) ||
    3.93 +             (type == GEQ && sum >= supply[n]) ||
    3.94 +             (type == LEQ && sum <= supply[n]);
    3.95 +    if (!b) return false;
    3.96    }
    3.97  
    3.98    return true;
    3.99 @@ -165,9 +177,10 @@
   3.100  // Check the feasibility of the given potentials (dual soluiton)
   3.101  // using the "Complementary Slackness" optimality condition
   3.102  template < typename GR, typename LM, typename UM,
   3.103 -           typename CM, typename FM, typename PM >
   3.104 +           typename CM, typename SM, typename FM, typename PM >
   3.105  bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
   3.106 -                     const CM& cost, const FM& flow, const PM& pi )
   3.107 +                     const CM& cost, const SM& supply, const FM& flow, 
   3.108 +                     const PM& pi )
   3.109  {
   3.110    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   3.111  
   3.112 @@ -179,6 +192,16 @@
   3.113            (red_cost > 0 && flow[e] == lower[e]) ||
   3.114            (red_cost < 0 && flow[e] == upper[e]);
   3.115    }
   3.116 +  
   3.117 +  for (NodeIt n(gr); opt && n != INVALID; ++n) {
   3.118 +    typename SM::Value sum = 0;
   3.119 +    for (OutArcIt e(gr, n); e != INVALID; ++e)
   3.120 +      sum += flow[e];
   3.121 +    for (InArcIt e(gr, n); e != INVALID; ++e)
   3.122 +      sum -= flow[e];
   3.123 +    opt = (sum == supply[n]) || (pi[n] == 0);
   3.124 +  }
   3.125 +  
   3.126    return opt;
   3.127  }
   3.128  
   3.129 @@ -190,14 +213,15 @@
   3.130                 const GR& gr, const LM& lower, const UM& upper,
   3.131                 const CM& cost, const SM& supply,
   3.132                 bool result, typename CM::Value total,
   3.133 -               const std::string &test_id = "" )
   3.134 +               const std::string &test_id = "",
   3.135 +               ProblemType type = EQ )
   3.136  {
   3.137    check(mcf_result == result, "Wrong result " + test_id);
   3.138    if (result) {
   3.139 -    check(checkFlow(gr, lower, upper, supply, mcf.flowMap()),
   3.140 +    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
   3.141            "The flow is not feasible " + test_id);
   3.142      check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
   3.143 -    check(checkPotential(gr, lower, upper, cost, mcf.flowMap(),
   3.144 +    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
   3.145                           mcf.potentialMap()),
   3.146            "Wrong potentials " + test_id);
   3.147    }
   3.148 @@ -226,7 +250,7 @@
   3.149    // Read the test digraph
   3.150    Digraph gr;
   3.151    Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
   3.152 -  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr);
   3.153 +  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
   3.154    ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
   3.155    Node v, w;
   3.156  
   3.157 @@ -239,6 +263,8 @@
   3.158      .nodeMap("sup1", s1)
   3.159      .nodeMap("sup2", s2)
   3.160      .nodeMap("sup3", s3)
   3.161 +    .nodeMap("sup4", s4)
   3.162 +    .nodeMap("sup5", s5)
   3.163      .node("source", v)
   3.164      .node("target", w)
   3.165      .run();
   3.166 @@ -247,6 +273,7 @@
   3.167    {
   3.168      NetworkSimplex<Digraph> mcf(gr);
   3.169  
   3.170 +    // Check the equality form
   3.171      mcf.upperMap(u).costMap(c);
   3.172      checkMcf(mcf, mcf.supplyMap(s1).run(),
   3.173               gr, l1, u, c, s1, true,  5240, "#A1");
   3.174 @@ -267,6 +294,28 @@
   3.175               gr, l1, cu, cc, s3, true,   0, "#A7");
   3.176      checkMcf(mcf, mcf.boundMaps(l2, u).run(),
   3.177               gr, l2, u, cc, s3, false,   0, "#A8");
   3.178 +
   3.179 +    // Check the GEQ form
   3.180 +    mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
   3.181 +    checkMcf(mcf, mcf.run(),
   3.182 +             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
   3.183 +    mcf.problemType(mcf.GEQ);
   3.184 +    checkMcf(mcf, mcf.lowerMap(l2).run(),
   3.185 +             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
   3.186 +    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
   3.187 +    checkMcf(mcf, mcf.run(),
   3.188 +             gr, l2, u, c, s5, false,    0, "#A11", GEQ);
   3.189 +
   3.190 +    // Check the LEQ form
   3.191 +    mcf.reset().problemType(mcf.LEQ);
   3.192 +    mcf.upperMap(u).costMap(c).supplyMap(s5);
   3.193 +    checkMcf(mcf, mcf.run(),
   3.194 +             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
   3.195 +    checkMcf(mcf, mcf.lowerMap(l2).run(),
   3.196 +             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
   3.197 +    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
   3.198 +    checkMcf(mcf, mcf.run(),
   3.199 +             gr, l2, u, c, s4, false,    0, "#A14", LEQ);
   3.200    }
   3.201  
   3.202    // B. Test NetworkSimplex with each pivot rule