1.1 --- a/lemon/circulation.h	Fri Apr 24 12:12:14 2009 +0100
     1.2 +++ b/lemon/circulation.h	Sat Apr 25 18:25:59 2009 +0200
     1.3 @@ -21,6 +21,7 @@
     1.4  
     1.5  #include <lemon/tolerance.h>
     1.6  #include <lemon/elevator.h>
     1.7 +#include <limits>
     1.8  
     1.9  ///\ingroup max_flow
    1.10  ///\file
    1.11 @@ -119,15 +120,15 @@
    1.12       at the nodes.
    1.13  
    1.14       The exact formulation of this problem is the following.
    1.15 -     Let \f$G=(V,A)\f$ be a digraph,
    1.16 -     \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$ denote the lower and
    1.17 -     upper bounds on the arcs, for which \f$0 \leq lower(uv) \leq upper(uv)\f$
    1.18 +     Let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$
    1.19 +     \f$upper: A\rightarrow\mathbf{R}\cup\{\infty\}\f$ denote the lower and
    1.20 +     upper bounds on the arcs, for which \f$lower(uv) \leq upper(uv)\f$
    1.21       holds for all \f$uv\in A\f$, and \f$sup: V\rightarrow\mathbf{R}\f$
    1.22       denotes the signed supply values of the nodes.
    1.23       If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
    1.24       supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
    1.25       \f$-sup(u)\f$ demand.
    1.26 -     A feasible circulation is an \f$f: A\rightarrow\mathbf{R}^+_0\f$
    1.27 +     A feasible circulation is an \f$f: A\rightarrow\mathbf{R}\f$
    1.28       solution of the following problem.
    1.29  
    1.30       \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu)
    1.31 @@ -151,6 +152,10 @@
    1.32       the direction of the arcs and taking the negative of the supply values
    1.33       (e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
    1.34  
    1.35 +     This algorithm either calculates a feasible circulation, or provides
    1.36 +     a \ref barrier() "barrier", which prooves that a feasible soultion
    1.37 +     cannot exist.
    1.38 +
    1.39       Note that this algorithm also provides a feasible solution for the
    1.40       \ref min_cost_flow "minimum cost flow problem".
    1.41  
    1.42 @@ -337,6 +342,13 @@
    1.43  
    1.44    private:
    1.45  
    1.46 +    bool checkBoundMaps() {
    1.47 +      for (ArcIt e(_g);e!=INVALID;++e) {
    1.48 +        if (_tol.less((*_up)[e], (*_lo)[e])) return false;
    1.49 +      }
    1.50 +      return true;
    1.51 +    }
    1.52 +
    1.53      void createStructures() {
    1.54        _node_num = _el = countNodes(_g);
    1.55  
    1.56 @@ -380,7 +392,7 @@
    1.57  
    1.58      /// Sets the upper bound (capacity) map.
    1.59      /// \return <tt>(*this)</tt>
    1.60 -    Circulation& upperMap(const LowerMap& map) {
    1.61 +    Circulation& upperMap(const UpperMap& map) {
    1.62        _up = &map;
    1.63        return *this;
    1.64      }
    1.65 @@ -467,6 +479,9 @@
    1.66      /// to the lower bound.
    1.67      void init()
    1.68      {
    1.69 +      LEMON_DEBUG(checkBoundMaps(),
    1.70 +        "Upper bounds must be greater or equal to the lower bounds");
    1.71 +
    1.72        createStructures();
    1.73  
    1.74        for(NodeIt n(_g);n!=INVALID;++n) {
    1.75 @@ -496,6 +511,9 @@
    1.76      /// to construct the initial solution.
    1.77      void greedyInit()
    1.78      {
    1.79 +      LEMON_DEBUG(checkBoundMaps(),
    1.80 +        "Upper bounds must be greater or equal to the lower bounds");
    1.81 +
    1.82        createStructures();
    1.83  
    1.84        for(NodeIt n(_g);n!=INVALID;++n) {
    1.85 @@ -503,11 +521,11 @@
    1.86        }
    1.87  
    1.88        for (ArcIt e(_g);e!=INVALID;++e) {
    1.89 -        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
    1.90 +        if (!_tol.less(-(*_excess)[_g.target(e)], (*_up)[e])) {
    1.91            _flow->set(e, (*_up)[e]);
    1.92            (*_excess)[_g.target(e)] += (*_up)[e];
    1.93            (*_excess)[_g.source(e)] -= (*_up)[e];
    1.94 -        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
    1.95 +        } else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) {
    1.96            _flow->set(e, (*_lo)[e]);
    1.97            (*_excess)[_g.target(e)] += (*_lo)[e];
    1.98            (*_excess)[_g.source(e)] -= (*_lo)[e];
    1.99 @@ -748,6 +766,9 @@
   1.100      bool checkBarrier() const
   1.101      {
   1.102        Flow delta=0;
   1.103 +      Flow inf_cap = std::numeric_limits<Flow>::has_infinity ?
   1.104 +        std::numeric_limits<Flow>::infinity() :
   1.105 +        std::numeric_limits<Flow>::max();
   1.106        for(NodeIt n(_g);n!=INVALID;++n)
   1.107          if(barrier(n))
   1.108            delta-=(*_supply)[n];
   1.109 @@ -755,7 +776,10 @@
   1.110          {
   1.111            Node s=_g.source(e);
   1.112            Node t=_g.target(e);
   1.113 -          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
   1.114 +          if(barrier(s)&&!barrier(t)) {
   1.115 +            if (_tol.less(inf_cap - (*_up)[e], delta)) return false;
   1.116 +            delta+=(*_up)[e];
   1.117 +          }
   1.118            else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
   1.119          }
   1.120        return _tol.negative(delta);