COIN-OR::LEMON - Graph Library

Ignore:
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • doc/groups.dox

    r658 r687  
    353353in a network with capacity constraints (lower and upper bounds)
    354354and arc costs.
    355 Formally, let \f$G=(V,A)\f$ be a digraph,
    356 \f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
     355Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
     356\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
    357357upper bounds for the flow values on the arcs, for which
    358 \f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
    359 \f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
    360 on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
     358\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
     359\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
     360on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
    361361signed supply values of the nodes.
    362362If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
    363363supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
    364364\f$-sup(u)\f$ demand.
    365 A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
     365A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
    366366of the following optimization problem.
    367367
     
    405405The dual solution of the minimum cost flow problem is represented by node
    406406potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
    407 An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
     407An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
    408408is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
    409409node potentials the following \e complementary \e slackness optimality
     
    414414   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
    415415   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
    416  - For all \f$u\in V\f$:
     416 - For all \f$u\in V\f$ nodes:
    417417   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
    418418     then \f$\pi(u)=0\f$.
    419419 
    420420Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
    421 \f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
     421\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
    422422\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
    423423
    424 All algorithms provide dual solution (node potentials) as well
     424All algorithms provide dual solution (node potentials) as well,
    425425if an optimal flow is found.
    426426
  • lemon/circulation.h

    r669 r688  
    6565    typedef SM SupplyMap;
    6666
    67     /// \brief The type of the flow values.
    68     typedef typename SupplyMap::Value Flow;
     67    /// \brief The type of the flow and supply values.
     68    typedef typename SupplyMap::Value Value;
    6969
    7070    /// \brief The type of the map that stores the flow values.
     
    7373    /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
    7474    /// concept.
    75     typedef typename Digraph::template ArcMap<Flow> FlowMap;
     75    typedef typename Digraph::template ArcMap<Value> FlowMap;
    7676
    7777    /// \brief Instantiates a FlowMap.
     
    105105    ///
    106106    /// The tolerance used by the algorithm to handle inexact computation.
    107     typedef lemon::Tolerance<Flow> Tolerance;
     107    typedef lemon::Tolerance<Value> Tolerance;
    108108
    109109  };
     
    188188    ///The type of the digraph the algorithm runs on.
    189189    typedef typename Traits::Digraph Digraph;
    190     ///The type of the flow values.
    191     typedef typename Traits::Flow Flow;
     190    ///The type of the flow and supply values.
     191    typedef typename Traits::Value Value;
    192192
    193193    ///The type of the lower bound map.
     
    222222    bool _local_level;
    223223
    224     typedef typename Digraph::template NodeMap<Flow> ExcessMap;
     224    typedef typename Digraph::template NodeMap<Value> ExcessMap;
    225225    ExcessMap* _excess;
    226226
     
    531531          (*_excess)[_g.source(e)] -= (*_lo)[e];
    532532        } else {
    533           Flow fc = -(*_excess)[_g.target(e)];
     533          Value fc = -(*_excess)[_g.target(e)];
    534534          _flow->set(e, fc);
    535535          (*_excess)[_g.target(e)] = 0;
     
    564564        int actlevel=(*_level)[act];
    565565        int mlevel=_node_num;
    566         Flow exc=(*_excess)[act];
     566        Value exc=(*_excess)[act];
    567567
    568568        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
    569569          Node v = _g.target(e);
    570           Flow fc=(*_up)[e]-(*_flow)[e];
     570          Value fc=(*_up)[e]-(*_flow)[e];
    571571          if(!_tol.positive(fc)) continue;
    572572          if((*_level)[v]<actlevel) {
     
    592592        for(InArcIt e(_g,act);e!=INVALID; ++e) {
    593593          Node v = _g.source(e);
    594           Flow fc=(*_flow)[e]-(*_lo)[e];
     594          Value fc=(*_flow)[e]-(*_lo)[e];
    595595          if(!_tol.positive(fc)) continue;
    596596          if((*_level)[v]<actlevel) {
     
    662662    ///@{
    663663
    664     /// \brief Returns the flow on the given arc.
    665     ///
    666     /// Returns the flow on the given arc.
     664    /// \brief Returns the flow value on the given arc.
     665    ///
     666    /// Returns the flow value on the given arc.
    667667    ///
    668668    /// \pre Either \ref run() or \ref init() must be called before
    669669    /// using this function.
    670     Flow flow(const Arc& arc) const {
     670    Value flow(const Arc& arc) const {
    671671      return (*_flow)[arc];
    672672    }
     
    751751      for(NodeIt n(_g);n!=INVALID;++n)
    752752        {
    753           Flow dif=-(*_supply)[n];
     753          Value dif=-(*_supply)[n];
    754754          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
    755755          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
     
    766766    bool checkBarrier() const
    767767    {
    768       Flow delta=0;
    769       Flow inf_cap = std::numeric_limits<Flow>::has_infinity ?
    770         std::numeric_limits<Flow>::infinity() :
    771         std::numeric_limits<Flow>::max();
     768      Value delta=0;
     769      Value inf_cap = std::numeric_limits<Value>::has_infinity ?
     770        std::numeric_limits<Value>::infinity() :
     771        std::numeric_limits<Value>::max();
    772772      for(NodeIt n(_g);n!=INVALID;++n)
    773773        if(barrier(n))
  • lemon/network_simplex.h

    r665 r690  
    3131#include <lemon/core.h>
    3232#include <lemon/math.h>
    33 #include <lemon/maps.h>
    34 #include <lemon/circulation.h>
    35 #include <lemon/adaptors.h>
    3633
    3734namespace lemon {
     
    5148  /// In general this class is the fastest implementation available
    5249  /// in LEMON for the minimum cost flow problem.
    53   /// Moreover it supports both direction of the supply/demand inequality
    54   /// constraints. For more information see \ref ProblemType.
     50  /// Moreover it supports both directions of the supply/demand inequality
     51  /// constraints. For more information see \ref SupplyType.
     52  ///
     53  /// Most of the parameters of the problem (except for the digraph)
     54  /// can be given using separate functions, and the algorithm can be
     55  /// executed using the \ref run() function. If some parameters are not
     56  /// specified, then default values will be used.
    5557  ///
    5658  /// \tparam GR The digraph type the algorithm runs on.
    57   /// \tparam F The value type used for flow amounts, capacity bounds
     59  /// \tparam V The value type used for flow amounts, capacity bounds
    5860  /// and supply values in the algorithm. By default it is \c int.
    5961  /// \tparam C The value type used for costs and potentials in the
    60   /// algorithm. By default it is the same as \c F.
     62  /// algorithm. By default it is the same as \c V.
    6163  ///
    6264  /// \warning Both value types must be signed and all input data must
     
    6668  /// implementations, from which the most efficient one is used
    6769  /// by default. For more information see \ref PivotRule.
    68   template <typename GR, typename F = int, typename C = F>
     70  template <typename GR, typename V = int, typename C = V>
    6971  class NetworkSimplex
    7072  {
    7173  public:
    7274
    73     /// The flow type of the algorithm
    74     typedef F Flow;
    75     /// The cost type of the algorithm
     75    /// The type of the flow amounts, capacity bounds and supply values
     76    typedef V Value;
     77    /// The type of the arc costs
    7678    typedef C Cost;
    77 #ifdef DOXYGEN
    78     /// The type of the flow map
    79     typedef GR::ArcMap<Flow> FlowMap;
    80     /// The type of the potential map
    81     typedef GR::NodeMap<Cost> PotentialMap;
    82 #else
    83     /// The type of the flow map
    84     typedef typename GR::template ArcMap<Flow> FlowMap;
    85     /// The type of the potential map
    86     typedef typename GR::template NodeMap<Cost> PotentialMap;
    87 #endif
    8879
    8980  public:
    9081
    91     /// \brief Enum type for selecting the pivot rule.
    92     ///
    93     /// Enum type for selecting the pivot rule for the \ref run()
    94     /// function.
    95     ///
    96     /// \ref NetworkSimplex provides five different pivot rule
    97     /// implementations that significantly affect the running time
    98     /// of the algorithm.
    99     /// By default \ref BLOCK_SEARCH "Block Search" is used, which
    100     /// proved to be the most efficient and the most robust on various
    101     /// test inputs according to our benchmark tests.
    102     /// However another pivot rule can be selected using the \ref run()
    103     /// function with the proper parameter.
    104     enum PivotRule {
    105 
    106       /// The First Eligible pivot rule.
    107       /// The next eligible arc is selected in a wraparound fashion
    108       /// in every iteration.
    109       FIRST_ELIGIBLE,
    110 
    111       /// The Best Eligible pivot rule.
    112       /// The best eligible arc is selected in every iteration.
    113       BEST_ELIGIBLE,
    114 
    115       /// The Block Search pivot rule.
    116       /// A specified number of arcs are examined in every iteration
    117       /// in a wraparound fashion and the best eligible arc is selected
    118       /// from this block.
    119       BLOCK_SEARCH,
    120 
    121       /// The Candidate List pivot rule.
    122       /// In a major iteration a candidate list is built from eligible arcs
    123       /// in a wraparound fashion and in the following minor iterations
    124       /// the best eligible arc is selected from this list.
    125       CANDIDATE_LIST,
    126 
    127       /// The Altering Candidate List pivot rule.
    128       /// It is a modified version of the Candidate List method.
    129       /// It keeps only the several best eligible arcs from the former
    130       /// candidate list and extends this list in every iteration.
    131       ALTERING_LIST
     82    /// \brief Problem type constants for the \c run() function.
     83    ///
     84    /// Enum type containing the problem type constants that can be
     85    /// returned by the \ref run() function of the algorithm.
     86    enum ProblemType {
     87      /// The problem has no feasible solution (flow).
     88      INFEASIBLE,
     89      /// The problem has optimal solution (i.e. it is feasible and
     90      /// bounded), and the algorithm has found optimal flow and node
     91      /// potentials (primal and dual solutions).
     92      OPTIMAL,
     93      /// The objective function of the problem is unbounded, i.e.
     94      /// there is a directed cycle having negative total cost and
     95      /// infinite upper bound.
     96      UNBOUNDED
    13297    };
    13398   
    134     /// \brief Enum type for selecting the problem type.
    135     ///
    136     /// Enum type for selecting the problem type, i.e. the direction of
    137     /// the inequalities in the supply/demand constraints of the
    138     /// \ref min_cost_flow "minimum cost flow problem".
    139     ///
    140     /// The default problem type is \c GEQ, since this form is supported
     99    /// \brief Constants for selecting the type of the supply constraints.
     100    ///
     101    /// Enum type containing constants for selecting the supply type,
     102    /// i.e. the direction of the inequalities in the supply/demand
     103    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
     104    ///
     105    /// The default supply type is \c GEQ, since this form is supported
    141106    /// by other minimum cost flow algorithms and the \ref Circulation
    142     /// algorithm as well.
    143     /// The \c LEQ problem type can be selected using the \ref problemType()
     107    /// algorithm, as well.
     108    /// The \c LEQ problem type can be selected using the \ref supplyType()
    144109    /// function.
    145110    ///
    146     /// Note that the equality form is a special case of both problem type.
    147     enum ProblemType {
    148 
    149       /// This option means that there are "<em>greater or equal</em>"
    150       /// constraints in the defintion, i.e. the exact formulation of the
    151       /// problem is the following.
     111    /// Note that the equality form is a special case of both supply types.
     112    enum SupplyType {
     113
     114      /// This option means that there are <em>"greater or equal"</em>
     115      /// supply/demand constraints in the definition, i.e. the exact
     116      /// formulation of the problem is the following.
    152117      /**
    153118          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
     
    165130      CARRY_SUPPLIES = GEQ,
    166131
    167       /// This option means that there are "<em>less or equal</em>"
    168       /// constraints in the defintion, i.e. the exact formulation of the
    169       /// problem is the following.
     132      /// This option means that there are <em>"less or equal"</em>
     133      /// supply/demand constraints in the definition, i.e. the exact
     134      /// formulation of the problem is the following.
    170135      /**
    171136          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
     
    183148      SATISFY_DEMANDS = LEQ
    184149    };
    185 
     150   
     151    /// \brief Constants for selecting the pivot rule.
     152    ///
     153    /// Enum type containing constants for selecting the pivot rule for
     154    /// the \ref run() function.
     155    ///
     156    /// \ref NetworkSimplex provides five different pivot rule
     157    /// implementations that significantly affect the running time
     158    /// of the algorithm.
     159    /// By default \ref BLOCK_SEARCH "Block Search" is used, which
     160    /// proved to be the most efficient and the most robust on various
     161    /// test inputs according to our benchmark tests.
     162    /// However another pivot rule can be selected using the \ref run()
     163    /// function with the proper parameter.
     164    enum PivotRule {
     165
     166      /// The First Eligible pivot rule.
     167      /// The next eligible arc is selected in a wraparound fashion
     168      /// in every iteration.
     169      FIRST_ELIGIBLE,
     170
     171      /// The Best Eligible pivot rule.
     172      /// The best eligible arc is selected in every iteration.
     173      BEST_ELIGIBLE,
     174
     175      /// The Block Search pivot rule.
     176      /// A specified number of arcs are examined in every iteration
     177      /// in a wraparound fashion and the best eligible arc is selected
     178      /// from this block.
     179      BLOCK_SEARCH,
     180
     181      /// The Candidate List pivot rule.
     182      /// In a major iteration a candidate list is built from eligible arcs
     183      /// in a wraparound fashion and in the following minor iterations
     184      /// the best eligible arc is selected from this list.
     185      CANDIDATE_LIST,
     186
     187      /// The Altering Candidate List pivot rule.
     188      /// It is a modified version of the Candidate List method.
     189      /// It keeps only the several best eligible arcs from the former
     190      /// candidate list and extends this list in every iteration.
     191      ALTERING_LIST
     192    };
     193   
    186194  private:
    187195
    188196    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    189 
    190     typedef typename GR::template ArcMap<Flow> FlowArcMap;
    191     typedef typename GR::template ArcMap<Cost> CostArcMap;
    192     typedef typename GR::template NodeMap<Flow> FlowNodeMap;
    193197
    194198    typedef std::vector<Arc> ArcVector;
     
    196200    typedef std::vector<int> IntVector;
    197201    typedef std::vector<bool> BoolVector;
    198     typedef std::vector<Flow> FlowVector;
     202    typedef std::vector<Value> ValueVector;
    199203    typedef std::vector<Cost> CostVector;
    200204
     
    214218
    215219    // Parameters of the problem
    216     FlowArcMap *_plower;
    217     FlowArcMap *_pupper;
    218     CostArcMap *_pcost;
    219     FlowNodeMap *_psupply;
    220     bool _pstsup;
    221     Node _psource, _ptarget;
    222     Flow _pstflow;
    223     ProblemType _ptype;
    224 
    225     // Result maps
    226     FlowMap *_flow_map;
    227     PotentialMap *_potential_map;
    228     bool _local_flow;
    229     bool _local_potential;
     220    bool _have_lower;
     221    SupplyType _stype;
     222    Value _sum_supply;
    230223
    231224    // Data structures for storing the digraph
    232225    IntNodeMap _node_id;
    233     ArcVector _arc_ref;
     226    IntArcMap _arc_id;
    234227    IntVector _source;
    235228    IntVector _target;
    236229
    237230    // Node and arc data
    238     FlowVector _cap;
     231    ValueVector _lower;
     232    ValueVector _upper;
     233    ValueVector _cap;
    239234    CostVector _cost;
    240     FlowVector _supply;
    241     FlowVector _flow;
     235    ValueVector _supply;
     236    ValueVector _flow;
    242237    CostVector _pi;
    243238
     
    258253    int first, second, right, last;
    259254    int stem, par_stem, new_stem;
    260     Flow delta;
     255    Value delta;
     256
     257  public:
     258 
     259    /// \brief Constant for infinite upper bounds (capacities).
     260    ///
     261    /// Constant for infinite upper bounds (capacities).
     262    /// It is \c std::numeric_limits<Value>::infinity() if available,
     263    /// \c std::numeric_limits<Value>::max() otherwise.
     264    const Value INF;
    261265
    262266  private:
     
    660664    /// \param graph The digraph the algorithm runs on.
    661665    NetworkSimplex(const GR& graph) :
    662       _graph(graph),
    663       _plower(NULL), _pupper(NULL), _pcost(NULL),
    664       _psupply(NULL), _pstsup(false), _ptype(GEQ),
    665       _flow_map(NULL), _potential_map(NULL),
    666       _local_flow(false), _local_potential(false),
    667       _node_id(graph)
     666      _graph(graph), _node_id(graph), _arc_id(graph),
     667      INF(std::numeric_limits<Value>::has_infinity ?
     668          std::numeric_limits<Value>::infinity() :
     669          std::numeric_limits<Value>::max())
    668670    {
    669       LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
    670                    std::numeric_limits<Flow>::is_signed,
    671         "The flow type of NetworkSimplex must be signed integer");
    672       LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
    673                    std::numeric_limits<Cost>::is_signed,
    674         "The cost type of NetworkSimplex must be signed integer");
    675     }
    676 
    677     /// Destructor.
    678     ~NetworkSimplex() {
    679       if (_local_flow) delete _flow_map;
    680       if (_local_potential) delete _potential_map;
     671      // Check the value types
     672      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
     673        "The flow type of NetworkSimplex must be signed");
     674      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
     675        "The cost type of NetworkSimplex must be signed");
     676       
     677      // Resize vectors
     678      _node_num = countNodes(_graph);
     679      _arc_num = countArcs(_graph);
     680      int all_node_num = _node_num + 1;
     681      int all_arc_num = _arc_num + _node_num;
     682
     683      _source.resize(all_arc_num);
     684      _target.resize(all_arc_num);
     685
     686      _lower.resize(all_arc_num);
     687      _upper.resize(all_arc_num);
     688      _cap.resize(all_arc_num);
     689      _cost.resize(all_arc_num);
     690      _supply.resize(all_node_num);
     691      _flow.resize(all_arc_num);
     692      _pi.resize(all_node_num);
     693
     694      _parent.resize(all_node_num);
     695      _pred.resize(all_node_num);
     696      _forward.resize(all_node_num);
     697      _thread.resize(all_node_num);
     698      _rev_thread.resize(all_node_num);
     699      _succ_num.resize(all_node_num);
     700      _last_succ.resize(all_node_num);
     701      _state.resize(all_arc_num);
     702
     703      // Copy the graph (store the arcs in a mixed order)
     704      int i = 0;
     705      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
     706        _node_id[n] = i;
     707      }
     708      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
     709      i = 0;
     710      for (ArcIt a(_graph); a != INVALID; ++a) {
     711        _arc_id[a] = i;
     712        _source[i] = _node_id[_graph.source(a)];
     713        _target[i] = _node_id[_graph.target(a)];
     714        if ((i += k) >= _arc_num) i = (i % k) + 1;
     715      }
     716     
     717      // Initialize maps
     718      for (int i = 0; i != _node_num; ++i) {
     719        _supply[i] = 0;
     720      }
     721      for (int i = 0; i != _arc_num; ++i) {
     722        _lower[i] = 0;
     723        _upper[i] = INF;
     724        _cost[i] = 1;
     725      }
     726      _have_lower = false;
     727      _stype = GEQ;
    681728    }
    682729
     
    690737    ///
    691738    /// This function sets the lower bounds on the arcs.
    692     /// If neither this function nor \ref boundMaps() is used before
    693     /// calling \ref run(), the lower bounds will be set to zero
    694     /// on all arcs.
     739    /// If it is not used before calling \ref run(), the lower bounds
     740    /// will be set to zero on all arcs.
    695741    ///
    696742    /// \param map An arc map storing the lower bounds.
    697     /// Its \c Value type must be convertible to the \c Flow type
     743    /// Its \c Value type must be convertible to the \c Value type
    698744    /// of the algorithm.
    699745    ///
    700746    /// \return <tt>(*this)</tt>
    701     template <typename LOWER>
    702     NetworkSimplex& lowerMap(const LOWER& map) {
    703       delete _plower;
    704       _plower = new FlowArcMap(_graph);
     747    template <typename LowerMap>
     748    NetworkSimplex& lowerMap(const LowerMap& map) {
     749      _have_lower = true;
    705750      for (ArcIt a(_graph); a != INVALID; ++a) {
    706         (*_plower)[a] = map[a];
     751        _lower[_arc_id[a]] = map[a];
    707752      }
    708753      return *this;
     
    712757    ///
    713758    /// This function sets the upper bounds (capacities) on the arcs.
    714     /// If none of the functions \ref upperMap(), \ref capacityMap()
    715     /// and \ref boundMaps() is used before calling \ref run(),
    716     /// the upper bounds (capacities) will be set to
    717     /// \c std::numeric_limits<Flow>::max() on all arcs.
     759    /// If it is not used before calling \ref run(), the upper bounds
     760    /// will be set to \ref INF on all arcs (i.e. the flow value will be
     761    /// unbounded from above on each arc).
    718762    ///
    719763    /// \param map An arc map storing the upper bounds.
    720     /// Its \c Value type must be convertible to the \c Flow type
     764    /// Its \c Value type must be convertible to the \c Value type
    721765    /// of the algorithm.
    722766    ///
    723767    /// \return <tt>(*this)</tt>
    724     template<typename UPPER>
    725     NetworkSimplex& upperMap(const UPPER& map) {
    726       delete _pupper;
    727       _pupper = new FlowArcMap(_graph);
     768    template<typename UpperMap>
     769    NetworkSimplex& upperMap(const UpperMap& map) {
    728770      for (ArcIt a(_graph); a != INVALID; ++a) {
    729         (*_pupper)[a] = map[a];
     771        _upper[_arc_id[a]] = map[a];
    730772      }
    731773      return *this;
    732     }
    733 
    734     /// \brief Set the upper bounds (capacities) on the arcs.
    735     ///
    736     /// This function sets the upper bounds (capacities) on the arcs.
    737     /// It is just an alias for \ref upperMap().
    738     ///
    739     /// \return <tt>(*this)</tt>
    740     template<typename CAP>
    741     NetworkSimplex& capacityMap(const CAP& map) {
    742       return upperMap(map);
    743     }
    744 
    745     /// \brief Set the lower and upper bounds on the arcs.
    746     ///
    747     /// This function sets the lower and upper bounds on the arcs.
    748     /// If neither this function nor \ref lowerMap() is used before
    749     /// calling \ref run(), the lower bounds will be set to zero
    750     /// on all arcs.
    751     /// If none of the functions \ref upperMap(), \ref capacityMap()
    752     /// and \ref boundMaps() is used before calling \ref run(),
    753     /// the upper bounds (capacities) will be set to
    754     /// \c std::numeric_limits<Flow>::max() on all arcs.
    755     ///
    756     /// \param lower An arc map storing the lower bounds.
    757     /// \param upper An arc map storing the upper bounds.
    758     ///
    759     /// The \c Value type of the maps must be convertible to the
    760     /// \c Flow type of the algorithm.
    761     ///
    762     /// \note This function is just a shortcut of calling \ref lowerMap()
    763     /// and \ref upperMap() separately.
    764     ///
    765     /// \return <tt>(*this)</tt>
    766     template <typename LOWER, typename UPPER>
    767     NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
    768       return lowerMap(lower).upperMap(upper);
    769774    }
    770775
     
    780785    ///
    781786    /// \return <tt>(*this)</tt>
    782     template<typename COST>
    783     NetworkSimplex& costMap(const COST& map) {
    784       delete _pcost;
    785       _pcost = new CostArcMap(_graph);
     787    template<typename CostMap>
     788    NetworkSimplex& costMap(const CostMap& map) {
    786789      for (ArcIt a(_graph); a != INVALID; ++a) {
    787         (*_pcost)[a] = map[a];
     790        _cost[_arc_id[a]] = map[a];
    788791      }
    789792      return *this;
     
    798801    ///
    799802    /// \param map A node map storing the supply values.
    800     /// Its \c Value type must be convertible to the \c Flow type
     803    /// Its \c Value type must be convertible to the \c Value type
    801804    /// of the algorithm.
    802805    ///
    803806    /// \return <tt>(*this)</tt>
    804     template<typename SUP>
    805     NetworkSimplex& supplyMap(const SUP& map) {
    806       delete _psupply;
    807       _pstsup = false;
    808       _psupply = new FlowNodeMap(_graph);
     807    template<typename SupplyMap>
     808    NetworkSimplex& supplyMap(const SupplyMap& map) {
    809809      for (NodeIt n(_graph); n != INVALID; ++n) {
    810         (*_psupply)[n] = map[n];
     810        _supply[_node_id[n]] = map[n];
    811811      }
    812812      return *this;
     
    821821    /// (It makes sense only if non-zero lower bounds are given.)
    822822    ///
     823    /// Using this function has the same effect as using \ref supplyMap()
     824    /// with such a map in which \c k is assigned to \c s, \c -k is
     825    /// assigned to \c t and all other nodes have zero supply value.
     826    ///
    823827    /// \param s The source node.
    824828    /// \param t The target node.
     
    827831    ///
    828832    /// \return <tt>(*this)</tt>
    829     NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
    830       delete _psupply;
    831       _psupply = NULL;
    832       _pstsup = true;
    833       _psource = s;
    834       _ptarget = t;
    835       _pstflow = k;
     833    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
     834      for (int i = 0; i != _node_num; ++i) {
     835        _supply[i] = 0;
     836      }
     837      _supply[_node_id[s]] =  k;
     838      _supply[_node_id[t]] = -k;
    836839      return *this;
    837840    }
    838841   
    839     /// \brief Set the problem type.
    840     ///
    841     /// This function sets the problem type for the algorithm.
    842     /// If it is not used before calling \ref run(), the \ref GEQ problem
     842    /// \brief Set the type of the supply constraints.
     843    ///
     844    /// This function sets the type of the supply/demand constraints.
     845    /// If it is not used before calling \ref run(), the \ref GEQ supply
    843846    /// type will be used.
    844847    ///
    845     /// For more information see \ref ProblemType.
     848    /// For more information see \ref SupplyType.
    846849    ///
    847850    /// \return <tt>(*this)</tt>
    848     NetworkSimplex& problemType(ProblemType problem_type) {
    849       _ptype = problem_type;
     851    NetworkSimplex& supplyType(SupplyType supply_type) {
     852      _stype = supply_type;
    850853      return *this;
    851854    }
    852855
    853     /// \brief Set the flow map.
    854     ///
    855     /// This function sets the flow map.
    856     /// If it is not used before calling \ref run(), an instance will
    857     /// be allocated automatically. The destructor deallocates this
    858     /// automatically allocated map, of course.
    859     ///
    860     /// \return <tt>(*this)</tt>
    861     NetworkSimplex& flowMap(FlowMap& map) {
    862       if (_local_flow) {
    863         delete _flow_map;
    864         _local_flow = false;
    865       }
    866       _flow_map = &map;
    867       return *this;
    868     }
    869 
    870     /// \brief Set the potential map.
    871     ///
    872     /// This function sets the potential map, which is used for storing
    873     /// the dual solution.
    874     /// If it is not used before calling \ref run(), an instance will
    875     /// be allocated automatically. The destructor deallocates this
    876     /// automatically allocated map, of course.
    877     ///
    878     /// \return <tt>(*this)</tt>
    879     NetworkSimplex& potentialMap(PotentialMap& map) {
    880       if (_local_potential) {
    881         delete _potential_map;
    882         _local_potential = false;
    883       }
    884       _potential_map = &map;
    885       return *this;
    886     }
    887    
    888856    /// @}
    889857
     
    897865    /// This function runs the algorithm.
    898866    /// The paramters can be specified using functions \ref lowerMap(),
    899     /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
    900     /// \ref costMap(), \ref supplyMap(), \ref stSupply(),
    901     /// \ref problemType(), \ref flowMap() and \ref potentialMap().
     867    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
     868    /// \ref supplyType().
    902869    /// For example,
    903870    /// \code
    904871    ///   NetworkSimplex<ListDigraph> ns(graph);
    905     ///   ns.boundMaps(lower, upper).costMap(cost)
     872    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
    906873    ///     .supplyMap(sup).run();
    907874    /// \endcode
     
    911878    /// \ref reset() is called, thus only the modified parameters
    912879    /// have to be set again. See \ref reset() for examples.
     880    /// However the underlying digraph must not be modified after this
     881    /// class have been constructed, since it copies and extends the graph.
    913882    ///
    914883    /// \param pivot_rule The pivot rule that will be used during the
    915884    /// algorithm. For more information see \ref PivotRule.
    916885    ///
    917     /// \return \c true if a feasible flow can be found.
    918     bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
    919       return init() && start(pivot_rule);
     886    /// \return \c INFEASIBLE if no feasible flow exists,
     887    /// \n \c OPTIMAL if the problem has optimal solution
     888    /// (i.e. it is feasible and bounded), and the algorithm has found
     889    /// optimal flow and node potentials (primal and dual solutions),
     890    /// \n \c UNBOUNDED if the objective function of the problem is
     891    /// unbounded, i.e. there is a directed cycle having negative total
     892    /// cost and infinite upper bound.
     893    ///
     894    /// \see ProblemType, PivotRule
     895    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
     896      if (!init()) return INFEASIBLE;
     897      return start(pivot_rule);
    920898    }
    921899
     
    924902    /// This function resets all the paramaters that have been given
    925903    /// before using functions \ref lowerMap(), \ref upperMap(),
    926     /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
    927     /// \ref supplyMap(), \ref stSupply(), \ref problemType(),
    928     /// \ref flowMap() and \ref potentialMap().
     904    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
    929905    ///
    930906    /// It is useful for multiple run() calls. If this function is not
    931907    /// used, all the parameters given before are kept for the next
    932908    /// \ref run() call.
     909    /// However the underlying digraph must not be modified after this
     910    /// class have been constructed, since it copies and extends the graph.
    933911    ///
    934912    /// For example,
     
    937915    ///
    938916    ///   // First run
    939     ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
     917    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
    940918    ///     .supplyMap(sup).run();
    941919    ///
     
    948926    ///   // (the lower bounds will be set to zero on all arcs)
    949927    ///   ns.reset();
    950     ///   ns.capacityMap(cap).costMap(cost)
     928    ///   ns.upperMap(capacity).costMap(cost)
    951929    ///     .supplyMap(sup).run();
    952930    /// \endcode
     
    954932    /// \return <tt>(*this)</tt>
    955933    NetworkSimplex& reset() {
    956       delete _plower;
    957       delete _pupper;
    958       delete _pcost;
    959       delete _psupply;
    960       _plower = NULL;
    961       _pupper = NULL;
    962       _pcost = NULL;
    963       _psupply = NULL;
    964       _pstsup = false;
    965       _ptype = GEQ;
    966       if (_local_flow) delete _flow_map;
    967       if (_local_potential) delete _potential_map;
    968       _flow_map = NULL;
    969       _potential_map = NULL;
    970       _local_flow = false;
    971       _local_potential = false;
    972 
     934      for (int i = 0; i != _node_num; ++i) {
     935        _supply[i] = 0;
     936      }
     937      for (int i = 0; i != _arc_num; ++i) {
     938        _lower[i] = 0;
     939        _upper[i] = INF;
     940        _cost[i] = 1;
     941      }
     942      _have_lower = false;
     943      _stype = GEQ;
    973944      return *this;
    974945    }
     
    986957    ///
    987958    /// This function returns the total cost of the found flow.
    988     /// The complexity of the function is O(e).
     959    /// Its complexity is O(e).
    989960    ///
    990961    /// \note The return type of the function can be specified as a
     
    998969    ///
    999970    /// \pre \ref run() must be called before using this function.
    1000     template <typename Num>
    1001     Num totalCost() const {
    1002       Num c = 0;
    1003       if (_pcost) {
    1004         for (ArcIt e(_graph); e != INVALID; ++e)
    1005           c += (*_flow_map)[e] * (*_pcost)[e];
    1006       } else {
    1007         for (ArcIt e(_graph); e != INVALID; ++e)
    1008           c += (*_flow_map)[e];
     971    template <typename Number>
     972    Number totalCost() const {
     973      Number c = 0;
     974      for (ArcIt a(_graph); a != INVALID; ++a) {
     975        int i = _arc_id[a];
     976        c += Number(_flow[i]) * Number(_cost[i]);
    1009977      }
    1010978      return c;
     
    1022990    ///
    1023991    /// \pre \ref run() must be called before using this function.
    1024     Flow flow(const Arc& a) const {
    1025       return (*_flow_map)[a];
    1026     }
    1027 
    1028     /// \brief Return a const reference to the flow map.
    1029     ///
    1030     /// This function returns a const reference to an arc map storing
    1031     /// the found flow.
     992    Value flow(const Arc& a) const {
     993      return _flow[_arc_id[a]];
     994    }
     995
     996    /// \brief Return the flow map (the primal solution).
     997    ///
     998    /// This function copies the flow value on each arc into the given
     999    /// map. The \c Value type of the algorithm must be convertible to
     1000    /// the \c Value type of the map.
    10321001    ///
    10331002    /// \pre \ref run() must be called before using this function.
    1034     const FlowMap& flowMap() const {
    1035       return *_flow_map;
     1003    template <typename FlowMap>
     1004    void flowMap(FlowMap &map) const {
     1005      for (ArcIt a(_graph); a != INVALID; ++a) {
     1006        map.set(a, _flow[_arc_id[a]]);
     1007      }
    10361008    }
    10371009
     
    10431015    /// \pre \ref run() must be called before using this function.
    10441016    Cost potential(const Node& n) const {
    1045       return (*_potential_map)[n];
    1046     }
    1047 
    1048     /// \brief Return a const reference to the potential map
    1049     /// (the dual solution).
    1050     ///
    1051     /// This function returns a const reference to a node map storing
    1052     /// the found potentials, which form the dual solution of the
    1053     /// \ref min_cost_flow "minimum cost flow" problem.
     1017      return _pi[_node_id[n]];
     1018    }
     1019
     1020    /// \brief Return the potential map (the dual solution).
     1021    ///
     1022    /// This function copies the potential (dual value) of each node
     1023    /// into the given map.
     1024    /// The \c Cost type of the algorithm must be convertible to the
     1025    /// \c Value type of the map.
    10541026    ///
    10551027    /// \pre \ref run() must be called before using this function.
    1056     const PotentialMap& potentialMap() const {
    1057       return *_potential_map;
     1028    template <typename PotentialMap>
     1029    void potentialMap(PotentialMap &map) const {
     1030      for (NodeIt n(_graph); n != INVALID; ++n) {
     1031        map.set(n, _pi[_node_id[n]]);
     1032      }
    10581033    }
    10591034
     
    10641039    // Initialize internal data structures
    10651040    bool init() {
    1066       // Initialize result maps
    1067       if (!_flow_map) {
    1068         _flow_map = new FlowMap(_graph);
    1069         _local_flow = true;
    1070       }
    1071       if (!_potential_map) {
    1072         _potential_map = new PotentialMap(_graph);
    1073         _local_potential = true;
    1074       }
    1075 
    1076       // Initialize vectors
    1077       _node_num = countNodes(_graph);
    1078       _arc_num = countArcs(_graph);
    1079       int all_node_num = _node_num + 1;
    1080       int all_arc_num = _arc_num + _node_num;
    10811041      if (_node_num == 0) return false;
    10821042
    1083       _arc_ref.resize(_arc_num);
    1084       _source.resize(all_arc_num);
    1085       _target.resize(all_arc_num);
    1086 
    1087       _cap.resize(all_arc_num);
    1088       _cost.resize(all_arc_num);
    1089       _supply.resize(all_node_num);
    1090       _flow.resize(all_arc_num);
    1091       _pi.resize(all_node_num);
    1092 
    1093       _parent.resize(all_node_num);
    1094       _pred.resize(all_node_num);
    1095       _forward.resize(all_node_num);
    1096       _thread.resize(all_node_num);
    1097       _rev_thread.resize(all_node_num);
    1098       _succ_num.resize(all_node_num);
    1099       _last_succ.resize(all_node_num);
    1100       _state.resize(all_arc_num);
    1101 
    1102       // Initialize node related data
    1103       bool valid_supply = true;
    1104       Flow sum_supply = 0;
    1105       if (!_pstsup && !_psupply) {
    1106         _pstsup = true;
    1107         _psource = _ptarget = NodeIt(_graph);
    1108         _pstflow = 0;
    1109       }
    1110       if (_psupply) {
    1111         int i = 0;
    1112         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
    1113           _node_id[n] = i;
    1114           _supply[i] = (*_psupply)[n];
    1115           sum_supply += _supply[i];
    1116         }
    1117         valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
    1118                        (_ptype == LEQ && sum_supply >= 0);
     1043      // Check the sum of supply values
     1044      _sum_supply = 0;
     1045      for (int i = 0; i != _node_num; ++i) {
     1046        _sum_supply += _supply[i];
     1047      }
     1048      if ( !((_stype == GEQ && _sum_supply <= 0) ||
     1049             (_stype == LEQ && _sum_supply >= 0)) ) return false;
     1050
     1051      // Remove non-zero lower bounds
     1052      if (_have_lower) {
     1053        for (int i = 0; i != _arc_num; ++i) {
     1054          Value c = _lower[i];
     1055          if (c >= 0) {
     1056            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
     1057          } else {
     1058            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
     1059          }
     1060          _supply[_source[i]] -= c;
     1061          _supply[_target[i]] += c;
     1062        }
    11191063      } else {
    1120         int i = 0;
    1121         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
    1122           _node_id[n] = i;
    1123           _supply[i] = 0;
    1124         }
    1125         _supply[_node_id[_psource]] =  _pstflow;
    1126         _supply[_node_id[_ptarget]] = -_pstflow;
    1127       }
    1128       if (!valid_supply) return false;
    1129 
    1130       // Infinite capacity value
    1131       Flow inf_cap =
    1132         std::numeric_limits<Flow>::has_infinity ?
    1133         std::numeric_limits<Flow>::infinity() :
    1134         std::numeric_limits<Flow>::max();
     1064        for (int i = 0; i != _arc_num; ++i) {
     1065          _cap[i] = _upper[i];
     1066        }
     1067      }
    11351068
    11361069      // Initialize artifical cost
    1137       Cost art_cost;
     1070      Cost ART_COST;
    11381071      if (std::numeric_limits<Cost>::is_exact) {
    1139         art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
     1072        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
    11401073      } else {
    1141         art_cost = std::numeric_limits<Cost>::min();
     1074        ART_COST = std::numeric_limits<Cost>::min();
    11421075        for (int i = 0; i != _arc_num; ++i) {
    1143           if (_cost[i] > art_cost) art_cost = _cost[i];
    1144         }
    1145         art_cost = (art_cost + 1) * _node_num;
    1146       }
    1147 
    1148       // Run Circulation to check if a feasible solution exists
    1149       typedef ConstMap<Arc, Flow> ConstArcMap;
    1150       ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
    1151       FlowNodeMap *csup = NULL;
    1152       bool local_csup = false;
    1153       if (_psupply) {
    1154         csup = _psupply;
    1155       } else {
    1156         csup = new FlowNodeMap(_graph, 0);
    1157         (*csup)[_psource] =  _pstflow;
    1158         (*csup)[_ptarget] = -_pstflow;
    1159         local_csup = true;
    1160       }
    1161       bool circ_result = false;
    1162       if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
    1163         // GEQ problem type
    1164         if (_plower) {
    1165           if (_pupper) {
    1166             Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
    1167               circ(_graph, *_plower, *_pupper, *csup);
    1168             circ_result = circ.run();
    1169           } else {
    1170             Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
    1171               circ(_graph, *_plower, inf_arc_map, *csup);
    1172             circ_result = circ.run();
    1173           }
    1174         } else {
    1175           if (_pupper) {
    1176             Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
    1177               circ(_graph, zero_arc_map, *_pupper, *csup);
    1178             circ_result = circ.run();
    1179           } else {
    1180             Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
    1181               circ(_graph, zero_arc_map, inf_arc_map, *csup);
    1182             circ_result = circ.run();
    1183           }
    1184         }
    1185       } else {
    1186         // LEQ problem type
    1187         typedef ReverseDigraph<const GR> RevGraph;
    1188         typedef NegMap<FlowNodeMap> NegNodeMap;
    1189         RevGraph rgraph(_graph);
    1190         NegNodeMap neg_csup(*csup);
    1191         if (_plower) {
    1192           if (_pupper) {
    1193             Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
    1194               circ(rgraph, *_plower, *_pupper, neg_csup);
    1195             circ_result = circ.run();
    1196           } else {
    1197             Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
    1198               circ(rgraph, *_plower, inf_arc_map, neg_csup);
    1199             circ_result = circ.run();
    1200           }
    1201         } else {
    1202           if (_pupper) {
    1203             Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
    1204               circ(rgraph, zero_arc_map, *_pupper, neg_csup);
    1205             circ_result = circ.run();
    1206           } else {
    1207             Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
    1208               circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
    1209             circ_result = circ.run();
    1210           }
    1211         }
    1212       }
    1213       if (local_csup) delete csup;
    1214       if (!circ_result) return false;
    1215 
     1076          if (_cost[i] > ART_COST) ART_COST = _cost[i];
     1077        }
     1078        ART_COST = (ART_COST + 1) * _node_num;
     1079      }
     1080
     1081      // Initialize arc maps
     1082      for (int i = 0; i != _arc_num; ++i) {
     1083        _flow[i] = 0;
     1084        _state[i] = STATE_LOWER;
     1085      }
     1086     
    12161087      // Set data for the artificial root node
    12171088      _root = _node_num;
     
    12201091      _thread[_root] = 0;
    12211092      _rev_thread[0] = _root;
    1222       _succ_num[_root] = all_node_num;
     1093      _succ_num[_root] = _node_num + 1;
    12231094      _last_succ[_root] = _root - 1;
    1224       _supply[_root] = -sum_supply;
    1225       if (sum_supply < 0) {
    1226         _pi[_root] = -art_cost;
    1227       } else {
    1228         _pi[_root] = art_cost;
    1229       }
    1230 
    1231       // Store the arcs in a mixed order
    1232       int k = std::max(int(std::sqrt(double(_arc_num))), 10);
    1233       int i = 0;
    1234       for (ArcIt e(_graph); e != INVALID; ++e) {
    1235         _arc_ref[i] = e;
    1236         if ((i += k) >= _arc_num) i = (i % k) + 1;
    1237       }
    1238 
    1239       // Initialize arc maps
    1240       if (_pupper && _pcost) {
    1241         for (int i = 0; i != _arc_num; ++i) {
    1242           Arc e = _arc_ref[i];
    1243           _source[i] = _node_id[_graph.source(e)];
    1244           _target[i] = _node_id[_graph.target(e)];
    1245           _cap[i] = (*_pupper)[e];
    1246           _cost[i] = (*_pcost)[e];
    1247           _flow[i] = 0;
    1248           _state[i] = STATE_LOWER;
    1249         }
    1250       } else {
    1251         for (int i = 0; i != _arc_num; ++i) {
    1252           Arc e = _arc_ref[i];
    1253           _source[i] = _node_id[_graph.source(e)];
    1254           _target[i] = _node_id[_graph.target(e)];
    1255           _flow[i] = 0;
    1256           _state[i] = STATE_LOWER;
    1257         }
    1258         if (_pupper) {
    1259           for (int i = 0; i != _arc_num; ++i)
    1260             _cap[i] = (*_pupper)[_arc_ref[i]];
    1261         } else {
    1262           for (int i = 0; i != _arc_num; ++i)
    1263             _cap[i] = inf_cap;
    1264         }
    1265         if (_pcost) {
    1266           for (int i = 0; i != _arc_num; ++i)
    1267             _cost[i] = (*_pcost)[_arc_ref[i]];
    1268         } else {
    1269           for (int i = 0; i != _arc_num; ++i)
    1270             _cost[i] = 1;
    1271         }
    1272       }
    1273      
    1274       // Remove non-zero lower bounds
    1275       if (_plower) {
    1276         for (int i = 0; i != _arc_num; ++i) {
    1277           Flow c = (*_plower)[_arc_ref[i]];
    1278           if (c != 0) {
    1279             _cap[i] -= c;
    1280             _supply[_source[i]] -= c;
    1281             _supply[_target[i]] += c;
    1282           }
    1283         }
    1284       }
     1095      _supply[_root] = -_sum_supply;
     1096      _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
    12851097
    12861098      // Add artificial arcs and initialize the spanning tree data structure
    12871099      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
     1100        _parent[u] = _root;
     1101        _pred[u] = e;
    12881102        _thread[u] = u + 1;
    12891103        _rev_thread[u + 1] = u;
    12901104        _succ_num[u] = 1;
    12911105        _last_succ[u] = u;
    1292         _parent[u] = _root;
    1293         _pred[u] = e;
    1294         _cost[e] = art_cost;
    1295         _cap[e] = inf_cap;
     1106        _cost[e] = ART_COST;
     1107        _cap[e] = INF;
    12961108        _state[e] = STATE_TREE;
    1297         if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
     1109        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
    12981110          _flow[e] = _supply[u];
    12991111          _forward[u] = true;
    1300           _pi[u] = -art_cost + _pi[_root];
     1112          _pi[u] = -ART_COST + _pi[_root];
    13011113        } else {
    13021114          _flow[e] = -_supply[u];
    13031115          _forward[u] = false;
    1304           _pi[u] = art_cost + _pi[_root];
     1116          _pi[u] = ART_COST + _pi[_root];
    13051117        }
    13061118      }
     
    13371149      delta = _cap[in_arc];
    13381150      int result = 0;
    1339       Flow d;
     1151      Value d;
    13401152      int e;
    13411153
     
    13431155      for (int u = first; u != join; u = _parent[u]) {
    13441156        e = _pred[u];
    1345         d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
     1157        d = _forward[u] ?
     1158          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
    13461159        if (d < delta) {
    13471160          delta = d;
     
    13531166      for (int u = second; u != join; u = _parent[u]) {
    13541167        e = _pred[u];
    1355         d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
     1168        d = _forward[u] ?
     1169          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
    13561170        if (d <= delta) {
    13571171          delta = d;
     
    13751189      // Augment along the cycle
    13761190      if (delta > 0) {
    1377         Flow val = _state[in_arc] * delta;
     1191        Value val = _state[in_arc] * delta;
    13781192        _flow[in_arc] += val;
    13791193        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
     
    15271341
    15281342    // Execute the algorithm
    1529     bool start(PivotRule pivot_rule) {
     1343    ProblemType start(PivotRule pivot_rule) {
    15301344      // Select the pivot rule implementation
    15311345      switch (pivot_rule) {
     
    15411355          return start<AlteringListPivotRule>();
    15421356      }
    1543       return false;
     1357      return INFEASIBLE; // avoid warning
    15441358    }
    15451359
    15461360    template <typename PivotRuleImpl>
    1547     bool start() {
     1361    ProblemType start() {
    15481362      PivotRuleImpl pivot(*this);
    15491363
     
    15521366        findJoinNode();
    15531367        bool change = findLeavingArc();
     1368        if (delta >= INF) return UNBOUNDED;
    15541369        changeFlow(change);
    15551370        if (change) {
     
    15581373        }
    15591374      }
    1560 
    1561       // Copy flow values to _flow_map
    1562       if (_plower) {
     1375     
     1376      // Check feasibility
     1377      if (_sum_supply < 0) {
     1378        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
     1379          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
     1380        }
     1381      }
     1382      else if (_sum_supply > 0) {
     1383        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
     1384          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
     1385        }
     1386      }
     1387      else {
     1388        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
     1389          if (_flow[e] != 0) return INFEASIBLE;
     1390        }
     1391      }
     1392
     1393      // Transform the solution and the supply map to the original form
     1394      if (_have_lower) {
    15631395        for (int i = 0; i != _arc_num; ++i) {
    1564           Arc e = _arc_ref[i];
    1565           _flow_map->set(e, (*_plower)[e] + _flow[i]);
    1566         }
    1567       } else {
    1568         for (int i = 0; i != _arc_num; ++i) {
    1569           _flow_map->set(_arc_ref[i], _flow[i]);
    1570         }
    1571       }
    1572       // Copy potential values to _potential_map
    1573       for (NodeIt n(_graph); n != INVALID; ++n) {
    1574         _potential_map->set(n, _pi[_node_id[n]]);
    1575       }
    1576 
    1577       return true;
     1396          Value c = _lower[i];
     1397          if (c != 0) {
     1398            _flow[i] += c;
     1399            _supply[_source[i]] += c;
     1400            _supply[_target[i]] -= c;
     1401          }
     1402        }
     1403      }
     1404
     1405      return OPTIMAL;
    15781406    }
    15791407
  • lemon/preflow.h

    r658 r688  
    4747
    4848    /// \brief The type of the flow values.
    49     typedef typename CapacityMap::Value Flow;
     49    typedef typename CapacityMap::Value Value;
    5050
    5151    /// \brief The type of the map that stores the flow values.
     
    5353    /// The type of the map that stores the flow values.
    5454    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
    55     typedef typename Digraph::template ArcMap<Flow> FlowMap;
     55    typedef typename Digraph::template ArcMap<Value> FlowMap;
    5656
    5757    /// \brief Instantiates a FlowMap.
     
    8585    ///
    8686    /// The tolerance used by the algorithm to handle inexact computation.
    87     typedef lemon::Tolerance<Flow> Tolerance;
     87    typedef lemon::Tolerance<Value> Tolerance;
    8888
    8989  };
     
    126126    typedef typename Traits::CapacityMap CapacityMap;
    127127    ///The type of the flow values.
    128     typedef typename Traits::Flow Flow;
     128    typedef typename Traits::Value Value;
    129129
    130130    ///The type of the flow map.
     
    152152    bool _local_level;
    153153
    154     typedef typename Digraph::template NodeMap<Flow> ExcessMap;
     154    typedef typename Digraph::template NodeMap<Value> ExcessMap;
    155155    ExcessMap* _excess;
    156156
     
    471471
    472472      for (NodeIt n(_graph); n != INVALID; ++n) {
    473         Flow excess = 0;
     473        Value excess = 0;
    474474        for (InArcIt e(_graph, n); e != INVALID; ++e) {
    475475          excess += (*_flow)[e];
     
    520520
    521521      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
    522         Flow rem = (*_capacity)[e] - (*_flow)[e];
     522        Value rem = (*_capacity)[e] - (*_flow)[e];
    523523        if (_tolerance.positive(rem)) {
    524524          Node u = _graph.target(e);
     
    532532      }
    533533      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
    534         Flow rem = (*_flow)[e];
     534        Value rem = (*_flow)[e];
    535535        if (_tolerance.positive(rem)) {
    536536          Node v = _graph.source(e);
     
    565565
    566566        while (num > 0 && n != INVALID) {
    567           Flow excess = (*_excess)[n];
     567          Value excess = (*_excess)[n];
    568568          int new_level = _level->maxLevel();
    569569
    570570          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
    571             Flow rem = (*_capacity)[e] - (*_flow)[e];
     571            Value rem = (*_capacity)[e] - (*_flow)[e];
    572572            if (!_tolerance.positive(rem)) continue;
    573573            Node v = _graph.target(e);
     
    592592
    593593          for (InArcIt e(_graph, n); e != INVALID; ++e) {
    594             Flow rem = (*_flow)[e];
     594            Value rem = (*_flow)[e];
    595595            if (!_tolerance.positive(rem)) continue;
    596596            Node v = _graph.source(e);
     
    638638        num = _node_num * 20;
    639639        while (num > 0 && n != INVALID) {
    640           Flow excess = (*_excess)[n];
     640          Value excess = (*_excess)[n];
    641641          int new_level = _level->maxLevel();
    642642
    643643          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
    644             Flow rem = (*_capacity)[e] - (*_flow)[e];
     644            Value rem = (*_capacity)[e] - (*_flow)[e];
    645645            if (!_tolerance.positive(rem)) continue;
    646646            Node v = _graph.target(e);
     
    665665
    666666          for (InArcIt e(_graph, n); e != INVALID; ++e) {
    667             Flow rem = (*_flow)[e];
     667            Value rem = (*_flow)[e];
    668668            if (!_tolerance.positive(rem)) continue;
    669669            Node v = _graph.source(e);
     
    779779      Node n;
    780780      while ((n = _level->highestActive()) != INVALID) {
    781         Flow excess = (*_excess)[n];
     781        Value excess = (*_excess)[n];
    782782        int level = _level->highestActiveLevel();
    783783        int new_level = _level->maxLevel();
    784784
    785785        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
    786           Flow rem = (*_capacity)[e] - (*_flow)[e];
     786          Value rem = (*_capacity)[e] - (*_flow)[e];
    787787          if (!_tolerance.positive(rem)) continue;
    788788          Node v = _graph.target(e);
     
    807807
    808808        for (InArcIt e(_graph, n); e != INVALID; ++e) {
    809           Flow rem = (*_flow)[e];
     809          Value rem = (*_flow)[e];
    810810          if (!_tolerance.positive(rem)) continue;
    811811          Node v = _graph.source(e);
     
    898898    /// \pre Either \ref run() or \ref init() must be called before
    899899    /// using this function.
    900     Flow flowValue() const {
     900    Value flowValue() const {
    901901      return (*_excess)[_target];
    902902    }
    903903
    904     /// \brief Returns the flow on the given arc.
    905     ///
    906     /// Returns the flow on the given arc. This method can
     904    /// \brief Returns the flow value on the given arc.
     905    ///
     906    /// Returns the flow value on the given arc. This method can
    907907    /// be called after the second phase of the algorithm.
    908908    ///
    909909    /// \pre Either \ref run() or \ref init() must be called before
    910910    /// using this function.
    911     Flow flow(const Arc& arc) const {
     911    Value flow(const Arc& arc) const {
    912912      return (*_flow)[arc];
    913913    }
  • test/min_cost_flow_test.cc

    r662 r689  
    1919#include <iostream>
    2020#include <fstream>
     21#include <limits>
    2122
    2223#include <lemon/list_graph.h>
     
    3435char test_lgf[] =
    3536  "@nodes\n"
    36   "label  sup1 sup2 sup3 sup4 sup5\n"
    37   "    1    20   27    0   20   30\n"
    38   "    2    -4    0    0   -8   -3\n"
    39   "    3     0    0    0    0    0\n"
    40   "    4     0    0    0    0    0\n"
    41   "    5     9    0    0    6   11\n"
    42   "    6    -6    0    0   -5   -6\n"
    43   "    7     0    0    0    0    0\n"
    44   "    8     0    0    0    0    3\n"
    45   "    9     3    0    0    0    0\n"
    46   "   10    -2    0    0   -7   -2\n"
    47   "   11     0    0    0  -10    0\n"
    48   "   12   -20  -27    0  -30  -20\n"
    49   "\n"
     37  "label  sup1 sup2 sup3 sup4 sup5 sup6\n"
     38  "    1    20   27    0   30   20   30\n"
     39  "    2    -4    0    0    0   -8   -3\n"
     40  "    3     0    0    0    0    0    0\n"
     41  "    4     0    0    0    0    0    0\n"
     42  "    5     9    0    0    0    6   11\n"
     43  "    6    -6    0    0    0   -5   -6\n"
     44  "    7     0    0    0    0    0    0\n"
     45  "    8     0    0    0    0    0    3\n"
     46  "    9     3    0    0    0    0    0\n"
     47  "   10    -2    0    0    0   -7   -2\n"
     48  "   11     0    0    0    0  -10    0\n"
     49  "   12   -20  -27    0  -30  -30  -20\n"
     50  "\n"               
    5051  "@arcs\n"
    51   "       cost  cap low1 low2\n"
    52   " 1  2    70   11    0    8\n"
    53   " 1  3   150    3    0    1\n"
    54   " 1  4    80   15    0    2\n"
    55   " 2  8    80   12    0    0\n"
    56   " 3  5   140    5    0    3\n"
    57   " 4  6    60   10    0    1\n"
    58   " 4  7    80    2    0    0\n"
    59   " 4  8   110    3    0    0\n"
    60   " 5  7    60   14    0    0\n"
    61   " 5 11   120   12    0    0\n"
    62   " 6  3     0    3    0    0\n"
    63   " 6  9   140    4    0    0\n"
    64   " 6 10    90    8    0    0\n"
    65   " 7  1    30    5    0    0\n"
    66   " 8 12    60   16    0    4\n"
    67   " 9 12    50    6    0    0\n"
    68   "10 12    70   13    0    5\n"
    69   "10  2   100    7    0    0\n"
    70   "10  7    60   10    0    0\n"
    71   "11 10    20   14    0    6\n"
    72   "12 11    30   10    0    0\n"
     52  "       cost  cap low1 low2 low3\n"
     53  " 1  2    70   11    0    8    8\n"
     54  " 1  3   150    3    0    1    0\n"
     55  " 1  4    80   15    0    2    2\n"
     56  " 2  8    80   12    0    0    0\n"
     57  " 3  5   140    5    0    3    1\n"
     58  " 4  6    60   10    0    1    0\n"
     59  " 4  7    80    2    0    0    0\n"
     60  " 4  8   110    3    0    0    0\n"
     61  " 5  7    60   14    0    0    0\n"
     62  " 5 11   120   12    0    0    0\n"
     63  " 6  3     0    3    0    0    0\n"
     64  " 6  9   140    4    0    0    0\n"
     65  " 6 10    90    8    0    0    0\n"
     66  " 7  1    30    5    0    0   -5\n"
     67  " 8 12    60   16    0    4    3\n"
     68  " 9 12    50    6    0    0    0\n"
     69  "10 12    70   13    0    5    2\n"
     70  "10  2   100    7    0    0    0\n"
     71  "10  7    60   10    0    0   -3\n"
     72  "11 10    20   14    0    6  -20\n"
     73  "12 11    30   10    0    0  -10\n"
    7374  "\n"
    7475  "@attributes\n"
     
    7778
    7879
    79 enum ProblemType {
     80enum SupplyType {
    8081  EQ,
    8182  GEQ,
     
    8485
    8586// Check the interface of an MCF algorithm
    86 template <typename GR, typename Flow, typename Cost>
     87template <typename GR, typename Value, typename Cost>
    8788class McfClassConcept
    8889{
     
    9596
    9697      MCF mcf(g);
     98      const MCF& const_mcf = mcf;
    9799
    98100      b = mcf.reset()
    99101             .lowerMap(lower)
    100102             .upperMap(upper)
    101              .capacityMap(upper)
    102              .boundMaps(lower, upper)
    103103             .costMap(cost)
    104104             .supplyMap(sup)
    105105             .stSupply(n, n, k)
    106              .flowMap(flow)
    107              .potentialMap(pot)
    108106             .run();
    109      
    110       const MCF& const_mcf = mcf;
    111 
    112       const typename MCF::FlowMap &fm = const_mcf.flowMap();
    113       const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
    114 
    115       v = const_mcf.totalCost();
    116       double x = const_mcf.template totalCost<double>();
     107
     108      c = const_mcf.totalCost();
     109      x = const_mcf.template totalCost<double>();
    117110      v = const_mcf.flow(a);
    118       v = const_mcf.potential(n);
    119 
    120       ignore_unused_variable_warning(fm);
    121       ignore_unused_variable_warning(pm);
    122       ignore_unused_variable_warning(x);
     111      c = const_mcf.potential(n);
     112      const_mcf.flowMap(fm);
     113      const_mcf.potentialMap(pm);
    123114    }
    124115
    125116    typedef typename GR::Node Node;
    126117    typedef typename GR::Arc Arc;
    127     typedef concepts::ReadMap<Node, Flow> NM;
    128     typedef concepts::ReadMap<Arc, Flow> FAM;
     118    typedef concepts::ReadMap<Node, Value> NM;
     119    typedef concepts::ReadMap<Arc, Value> VAM;
    129120    typedef concepts::ReadMap<Arc, Cost> CAM;
     121    typedef concepts::WriteMap<Arc, Value> FlowMap;
     122    typedef concepts::WriteMap<Node, Cost> PotMap;
    130123
    131124    const GR &g;
    132     const FAM &lower;
    133     const FAM &upper;
     125    const VAM &lower;
     126    const VAM &upper;
    134127    const CAM &cost;
    135128    const NM &sup;
    136129    const Node &n;
    137130    const Arc &a;
    138     const Flow &k;
    139     Flow v;
     131    const Value &k;
     132    FlowMap fm;
     133    PotMap pm;
    140134    bool b;
    141 
    142     typename MCF::FlowMap &flow;
    143     typename MCF::PotentialMap &pot;
     135    double x;
     136    typename MCF::Value v;
     137    typename MCF::Cost c;
    144138  };
    145139
     
    152146bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
    153147                const SM& supply, const FM& flow,
    154                 ProblemType type = EQ )
     148                SupplyType type = EQ )
    155149{
    156150  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
     
    209203template < typename MCF, typename GR,
    210204           typename LM, typename UM,
    211            typename CM, typename SM >
    212 void checkMcf( const MCF& mcf, bool mcf_result,
     205           typename CM, typename SM,
     206           typename PT >
     207void checkMcf( const MCF& mcf, PT mcf_result,
    213208               const GR& gr, const LM& lower, const UM& upper,
    214209               const CM& cost, const SM& supply,
    215                bool result, typename CM::Value total,
     210               PT result, bool optimal, typename CM::Value total,
    216211               const std::string &test_id = "",
    217                ProblemType type = EQ )
     212               SupplyType type = EQ )
    218213{
    219214  check(mcf_result == result, "Wrong result " + test_id);
    220   if (result) {
    221     check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
     215  if (optimal) {
     216    typename GR::template ArcMap<typename SM::Value> flow(gr);
     217    typename GR::template NodeMap<typename CM::Value> pi(gr);
     218    mcf.flowMap(flow);
     219    mcf.potentialMap(pi);
     220    check(checkFlow(gr, lower, upper, supply, flow, type),
    222221          "The flow is not feasible " + test_id);
    223222    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
    224     check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
    225                          mcf.potentialMap()),
     223    check(checkPotential(gr, lower, upper, cost, supply, flow, pi),
    226224          "Wrong potentials " + test_id);
    227225  }
     
    232230  // Check the interfaces
    233231  {
    234     typedef int Flow;
    235     typedef int Cost;
    236232    typedef concepts::Digraph GR;
    237     checkConcept< McfClassConcept<GR, Flow, Cost>,
    238                   NetworkSimplex<GR, Flow, Cost> >();
     233    checkConcept< McfClassConcept<GR, int, int>,
     234                  NetworkSimplex<GR> >();
     235    checkConcept< McfClassConcept<GR, double, double>,
     236                  NetworkSimplex<GR, double> >();
     237    checkConcept< McfClassConcept<GR, int, double>,
     238                  NetworkSimplex<GR, int, double> >();
    239239  }
    240240
     
    245245  // Read the test digraph
    246246  Digraph gr;
    247   Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
    248   Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
     247  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
     248  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
    249249  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
    250250  Node v, w;
     
    256256    .arcMap("low1", l1)
    257257    .arcMap("low2", l2)
     258    .arcMap("low3", l3)
    258259    .nodeMap("sup1", s1)
    259260    .nodeMap("sup2", s2)
     
    261262    .nodeMap("sup4", s4)
    262263    .nodeMap("sup5", s5)
     264    .nodeMap("sup6", s6)
    263265    .node("source", v)
    264266    .node("target", w)
    265267    .run();
     268 
     269  // Build a test digraph for testing negative costs
     270  Digraph ngr;
     271  Node n1 = ngr.addNode();
     272  Node n2 = ngr.addNode();
     273  Node n3 = ngr.addNode();
     274  Node n4 = ngr.addNode();
     275  Node n5 = ngr.addNode();
     276  Node n6 = ngr.addNode();
     277  Node n7 = ngr.addNode();
     278 
     279  Arc a1 = ngr.addArc(n1, n2);
     280  Arc a2 = ngr.addArc(n1, n3);
     281  Arc a3 = ngr.addArc(n2, n4);
     282  Arc a4 = ngr.addArc(n3, n4);
     283  Arc a5 = ngr.addArc(n3, n2);
     284  Arc a6 = ngr.addArc(n5, n3);
     285  Arc a7 = ngr.addArc(n5, n6);
     286  Arc a8 = ngr.addArc(n6, n7);
     287  Arc a9 = ngr.addArc(n7, n5);
     288 
     289  Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
     290  ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
     291  Digraph::NodeMap<int> ns(ngr, 0);
     292 
     293  nl2[a7] =  1000;
     294  nl2[a8] = -1000;
     295 
     296  ns[n1] =  100;
     297  ns[n4] = -100;
     298 
     299  nc[a1] =  100;
     300  nc[a2] =   30;
     301  nc[a3] =   20;
     302  nc[a4] =   80;
     303  nc[a5] =   50;
     304  nc[a6] =   10;
     305  nc[a7] =   80;
     306  nc[a8] =   30;
     307  nc[a9] = -120;
    266308
    267309  // A. Test NetworkSimplex with the default pivot rule
     
    272314    mcf.upperMap(u).costMap(c);
    273315    checkMcf(mcf, mcf.supplyMap(s1).run(),
    274              gr, l1, u, c, s1, true,  5240, "#A1");
     316             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
    275317    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
    276              gr, l1, u, c, s2, true,  7620, "#A2");
     318             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
    277319    mcf.lowerMap(l2);
    278320    checkMcf(mcf, mcf.supplyMap(s1).run(),
    279              gr, l2, u, c, s1, true,  5970, "#A3");
     321             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
    280322    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
    281              gr, l2, u, c, s2, true,  8010, "#A4");
     323             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
    282324    mcf.reset();
    283325    checkMcf(mcf, mcf.supplyMap(s1).run(),
    284              gr, l1, cu, cc, s1, true,  74, "#A5");
     326             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
    285327    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
    286              gr, l2, cu, cc, s2, true,  94, "#A6");
     328             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
    287329    mcf.reset();
    288330    checkMcf(mcf, mcf.run(),
    289              gr, l1, cu, cc, s3, true,   0, "#A7");
    290     checkMcf(mcf, mcf.boundMaps(l2, u).run(),
    291              gr, l2, u, cc, s3, false,   0, "#A8");
     331             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
     332    checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
     333             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
     334    mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
     335    checkMcf(mcf, mcf.run(),
     336             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
    292337
    293338    // Check the GEQ form
    294     mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
    295     checkMcf(mcf, mcf.run(),
    296              gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
    297     mcf.problemType(mcf.GEQ);
     339    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
     340    checkMcf(mcf, mcf.run(),
     341             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
     342    mcf.supplyType(mcf.GEQ);
    298343    checkMcf(mcf, mcf.lowerMap(l2).run(),
    299              gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
    300     mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
    301     checkMcf(mcf, mcf.run(),
    302              gr, l2, u, c, s5, false,    0, "#A11", GEQ);
     344             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
     345    mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
     346    checkMcf(mcf, mcf.run(),
     347             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
    303348
    304349    // Check the LEQ form
    305     mcf.reset().problemType(mcf.LEQ);
    306     mcf.upperMap(u).costMap(c).supplyMap(s5);
    307     checkMcf(mcf, mcf.run(),
    308              gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
     350    mcf.reset().supplyType(mcf.LEQ);
     351    mcf.upperMap(u).costMap(c).supplyMap(s6);
     352    checkMcf(mcf, mcf.run(),
     353             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
    309354    checkMcf(mcf, mcf.lowerMap(l2).run(),
    310              gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
    311     mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
    312     checkMcf(mcf, mcf.run(),
    313              gr, l2, u, c, s4, false,    0, "#A14", LEQ);
     355             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
     356    mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
     357    checkMcf(mcf, mcf.run(),
     358             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
     359
     360    // Check negative costs
     361    NetworkSimplex<Digraph> nmcf(ngr);
     362    nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
     363    checkMcf(nmcf, nmcf.run(),
     364      ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
     365    checkMcf(nmcf, nmcf.upperMap(nu2).run(),
     366      ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
     367    nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
     368    checkMcf(nmcf, nmcf.run(),
     369      ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
    314370  }
    315371
     
    317373  {
    318374    NetworkSimplex<Digraph> mcf(gr);
    319     mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
     375    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
    320376
    321377    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
    322              gr, l2, u, c, s1, true,  5970, "#B1");
     378             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
    323379    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
    324              gr, l2, u, c, s1, true,  5970, "#B2");
     380             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
    325381    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
    326              gr, l2, u, c, s1, true,  5970, "#B3");
     382             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
    327383    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
    328              gr, l2, u, c, s1, true,  5970, "#B4");
     384             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
    329385    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
    330              gr, l2, u, c, s1, true,  5970, "#B5");
     386             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
    331387  }
    332388
  • tools/dimacs-solver.cc

    r674 r691  
    120120  ti.restart();
    121121  NetworkSimplex<Digraph, Value> ns(g);
    122   ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
    123   if (sum_sup > 0) ns.problemType(ns.LEQ);
     122  ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup);
     123  if (sum_sup > 0) ns.supplyType(ns.LEQ);
    124124  if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';
    125125  ti.restart();
Note: See TracChangeset for help on using the changeset viewer.