COIN-OR::LEMON - Graph Library

About

Downloads

Documentation

Development

Issue Tracker

Browse Tickets

Custom Query

Roadmap

Browse Source

Local Info

Mailing Lists

Search

Changes in / [639:72ac25ad276e:644:8d289c89d43e] in lemon-main


Ignore:
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • doc/groups.dox

    r611 r640  
    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

    r622 r641  
    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

    r618 r643  
    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

    r611 r641  
    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

    r615 r642  
    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

    r627 r644  
    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.