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 / [644:8d289c89d43e:639:72ac25ad276e] in lemon-1.2


Ignore:
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • doc/groups.dox

    r640 r611  
    353353in a network with capacity constraints (lower and upper bounds)
    354354and arc costs.
    355 Formally, 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
     355Formally, let \f$G=(V,A)\f$ be a digraph,
     356\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
    357357upper bounds for the flow values on the arcs, for which
    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
    360 on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
     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
     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}\f$ solution
     365A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\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}\f$ feasible solution of the problem
     407An \f$f: A\rightarrow\mathbf{Z}^+_0\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$ nodes:
     416 - For all \f$u\in V\f$:
    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 potential function \f$\pi\f$, i.e.
     421\f$uv\in A\f$ with respect to the node potentials \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

    r641 r622  
    6565    typedef SM SupplyMap;
    6666
    67     /// \brief The type of the flow and supply values.
    68     typedef typename SupplyMap::Value Value;
     67    /// \brief The type of the flow values.
     68    typedef typename SupplyMap::Value Flow;
    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<Value> FlowMap;
     75    typedef typename Digraph::template ArcMap<Flow> FlowMap;
    7676
    7777    /// \brief Instantiates a FlowMap.
    … …  
    105105    ///
    106106    /// The tolerance used by the algorithm to handle inexact computation.
    107     typedef lemon::Tolerance<Value> Tolerance;
     107    typedef lemon::Tolerance<Flow> 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 and supply values.
    191     typedef typename Traits::Value Value;
     190    ///The type of the flow values.
     191    typedef typename Traits::Flow Flow;
    192192
    193193    ///The type of the lower bound map.
    … …  
    222222    bool _local_level;
    223223
    224     typedef typename Digraph::template NodeMap<Value> ExcessMap;
     224    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
    225225    ExcessMap* _excess;
    226226
    … …  
    531531          (*_excess)[_g.source(e)] -= (*_lo)[e];
    532532        } else {
    533           Value fc = -(*_excess)[_g.target(e)];
     533          Flow 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         Value exc=(*_excess)[act];
     566        Flow exc=(*_excess)[act];
    567567
    568568        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
    569569          Node v = _g.target(e);
    570           Value fc=(*_up)[e]-(*_flow)[e];
     570          Flow 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           Value fc=(*_flow)[e]-(*_lo)[e];
     594          Flow fc=(*_flow)[e]-(*_lo)[e];
    595595          if(!_tol.positive(fc)) continue;
    596596          if((*_level)[v]<actlevel) {
    … …  
    662662    ///@{
    663663
    664     /// \brief Returns the flow value on the given arc.
    665     ///
    666     /// Returns the flow value on the given arc.
     664    /// \brief Returns the flow on the given arc.
     665    ///
     666    /// Returns the flow on the given arc.
    667667    ///
    668668    /// \pre Either \ref run() or \ref init() must be called before
    669669    /// using this function.
    670     Value flow(const Arc& arc) const {
     670    Flow flow(const Arc& arc) const {
    671671      return (*_flow)[arc];
    672672    }
    … …  
    751751      for(NodeIt n(_g);n!=INVALID;++n)
    752752        {
    753           Value dif=-(*_supply)[n];
     753          Flow 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       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();
     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();
    772772      for(NodeIt n(_g);n!=INVALID;++n)
    773773        if(barrier(n))

  • lemon/network_simplex.h

    r643 r618  
    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>
    3336
    3437namespace lemon {
    … …  
    4851  /// In general this class is the fastest implementation available
    4952  /// in LEMON for the minimum cost flow problem.
    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.
     53  /// Moreover it supports both direction of the supply/demand inequality
     54  /// constraints. For more information see \ref ProblemType.
    5755  ///
    5856  /// \tparam GR The digraph type the algorithm runs on.
    59   /// \tparam V The value type used for flow amounts, capacity bounds
     57  /// \tparam F The value type used for flow amounts, capacity bounds
    6058  /// and supply values in the algorithm. By default it is \c int.
    6159  /// \tparam C The value type used for costs and potentials in the
    62   /// algorithm. By default it is the same as \c V.
     60  /// algorithm. By default it is the same as \c F.
    6361  ///
    6462  /// \warning Both value types must be signed and all input data must
    … …  
    6866  /// implementations, from which the most efficient one is used
    6967  /// by default. For more information see \ref PivotRule.
    70   template <typename GR, typename V = int, typename C = V>
     68  template <typename GR, typename F = int, typename C = F>
    7169  class NetworkSimplex
    7270  {
    7371  public:
    7472
    75     /// The type of the flow amounts, capacity bounds and supply values
    76     typedef V Value;
    77     /// The type of the arc costs
     73    /// The flow type of the algorithm
     74    typedef F Flow;
     75    /// The cost type of the algorithm
    7876    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
    7988
    8089  public:
    8190
    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
     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
    97132    };
    98133   
    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
     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
    106141    /// by other minimum cost flow algorithms and the \ref Circulation
    107     /// algorithm, as well.
    108     /// The \c LEQ problem type can be selected using the \ref supplyType()
     142    /// algorithm as well.
     143    /// The \c LEQ problem type can be selected using the \ref problemType()
    109144    /// function.
    110145    ///
    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.
     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.
    117152      /**
    118153          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    … …  
    130165      CARRY_SUPPLIES = GEQ,
    131166
    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.
     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.
    135170      /**
    136171          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    … …  
    148183      SATISFY_DEMANDS = LEQ
    149184    };
    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    
     185
    194186  private:
    195187
    196188    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;
    197193
    198194    typedef std::vector<Arc> ArcVector;
    … …  
    200196    typedef std::vector<int> IntVector;
    201197    typedef std::vector<bool> BoolVector;
    202     typedef std::vector<Value> ValueVector;
     198    typedef std::vector<Flow> FlowVector;
    203199    typedef std::vector<Cost> CostVector;
    204200
    … …  
    218214
    219215    // Parameters of the problem
    220     bool _have_lower;
    221     SupplyType _stype;
    222     Value _sum_supply;
     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;
    223230
    224231    // Data structures for storing the digraph
    225232    IntNodeMap _node_id;
    226     IntArcMap _arc_id;
     233    ArcVector _arc_ref;
    227234    IntVector _source;
    228235    IntVector _target;
    229236
    230237    // Node and arc data
    231     ValueVector _lower;
    232     ValueVector _upper;
    233     ValueVector _cap;
     238    FlowVector _cap;
    234239    CostVector _cost;
    235     ValueVector _supply;
    236     ValueVector _flow;
     240    FlowVector _supply;
     241    FlowVector _flow;
    237242    CostVector _pi;
    238243
    … …  
    253258    int first, second, right, last;
    254259    int stem, par_stem, new_stem;
    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;
     260    Flow delta;
    265261
    266262  private:
    … …  
    664660    /// \param graph The digraph the algorithm runs on.
    665661    NetworkSimplex(const GR& 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())
     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)
    670668    {
    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
     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;
     681    }
     682
     683    /// \name Parameters
     684    /// The parameters of the algorithm can be specified using these
     685    /// functions.
     686
     687    /// @{
     688
     689    /// \brief Set the lower bounds on the arcs.
     690    ///
     691    /// 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.
     695    ///
     696    /// \param map An arc map storing the lower bounds.
     697    /// Its \c Value type must be convertible to the \c Flow type
     698    /// of the algorithm.
     699    ///
     700    /// \return <tt>(*this)</tt>
     701    template <typename LOWER>
     702    NetworkSimplex& lowerMap(const LOWER& map) {
     703      delete _plower;
     704      _plower = new FlowArcMap(_graph);
     705      for (ArcIt a(_graph); a != INVALID; ++a) {
     706        (*_plower)[a] = map[a];
     707      }
     708      return *this;
     709    }
     710
     711    /// \brief Set the upper bounds (capacities) on the arcs.
     712    ///
     713    /// 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.
     718    ///
     719    /// \param map An arc map storing the upper bounds.
     720    /// Its \c Value type must be convertible to the \c Flow type
     721    /// of the algorithm.
     722    ///
     723    /// \return <tt>(*this)</tt>
     724    template<typename UPPER>
     725    NetworkSimplex& upperMap(const UPPER& map) {
     726      delete _pupper;
     727      _pupper = new FlowArcMap(_graph);
     728      for (ArcIt a(_graph); a != INVALID; ++a) {
     729        (*_pupper)[a] = map[a];
     730      }
     731      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);
     769    }
     770
     771    /// \brief Set the costs of the arcs.
     772    ///
     773    /// This function sets the costs of the arcs.
     774    /// If it is not used before calling \ref run(), the costs
     775    /// will be set to \c 1 on all arcs.
     776    ///
     777    /// \param map An arc map storing the costs.
     778    /// Its \c Value type must be convertible to the \c Cost type
     779    /// of the algorithm.
     780    ///
     781    /// \return <tt>(*this)</tt>
     782    template<typename COST>
     783    NetworkSimplex& costMap(const COST& map) {
     784      delete _pcost;
     785      _pcost = new CostArcMap(_graph);
     786      for (ArcIt a(_graph); a != INVALID; ++a) {
     787        (*_pcost)[a] = map[a];
     788      }
     789      return *this;
     790    }
     791
     792    /// \brief Set the supply values of the nodes.
     793    ///
     794    /// This function sets the supply values of the nodes.
     795    /// If neither this function nor \ref stSupply() is used before
     796    /// calling \ref run(), the supply of each node will be set to zero.
     797    /// (It makes sense only if non-zero lower bounds are given.)
     798    ///
     799    /// \param map A node map storing the supply values.
     800    /// Its \c Value type must be convertible to the \c Flow type
     801    /// of the algorithm.
     802    ///
     803    /// \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);
     809      for (NodeIt n(_graph); n != INVALID; ++n) {
     810        (*_psupply)[n] = map[n];
     811      }
     812      return *this;
     813    }
     814
     815    /// \brief Set single source and target nodes and a supply value.
     816    ///
     817    /// This function sets a single source node and a single target node
     818    /// and the required flow value.
     819    /// If neither this function nor \ref supplyMap() is used before
     820    /// calling \ref run(), the supply of each node will be set to zero.
     821    /// (It makes sense only if non-zero lower bounds are given.)
     822    ///
     823    /// \param s The source node.
     824    /// \param t The target node.
     825    /// \param k The required amount of flow from node \c s to node \c t
     826    /// (i.e. the supply of \c s and the demand of \c t).
     827    ///
     828    /// \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;
     836      return *this;
     837    }
     838   
     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
     843    /// type will be used.
     844    ///
     845    /// For more information see \ref ProblemType.
     846    ///
     847    /// \return <tt>(*this)</tt>
     848    NetworkSimplex& problemType(ProblemType problem_type) {
     849      _ptype = problem_type;
     850      return *this;
     851    }
     852
     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   
     888    /// @}
     889
     890    /// \name Execution Control
     891    /// The algorithm can be executed using \ref run().
     892
     893    /// @{
     894
     895    /// \brief Run the algorithm.
     896    ///
     897    /// This function runs the algorithm.
     898    /// 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().
     902    /// For example,
     903    /// \code
     904    ///   NetworkSimplex<ListDigraph> ns(graph);
     905    ///   ns.boundMaps(lower, upper).costMap(cost)
     906    ///     .supplyMap(sup).run();
     907    /// \endcode
     908    ///
     909    /// This function can be called more than once. All the parameters
     910    /// that have been given are kept for the next call, unless
     911    /// \ref reset() is called, thus only the modified parameters
     912    /// have to be set again. See \ref reset() for examples.
     913    ///
     914    /// \param pivot_rule The pivot rule that will be used during the
     915    /// algorithm. For more information see \ref PivotRule.
     916    ///
     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);
     920    }
     921
     922    /// \brief Reset all the parameters that have been given before.
     923    ///
     924    /// This function resets all the paramaters that have been given
     925    /// 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().
     929    ///
     930    /// It is useful for multiple run() calls. If this function is not
     931    /// used, all the parameters given before are kept for the next
     932    /// \ref run() call.
     933    ///
     934    /// For example,
     935    /// \code
     936    ///   NetworkSimplex<ListDigraph> ns(graph);
     937    ///
     938    ///   // First run
     939    ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
     940    ///     .supplyMap(sup).run();
     941    ///
     942    ///   // Run again with modified cost map (reset() is not called,
     943    ///   // so only the cost map have to be set again)
     944    ///   cost[e] += 100;
     945    ///   ns.costMap(cost).run();
     946    ///
     947    ///   // Run again from scratch using reset()
     948    ///   // (the lower bounds will be set to zero on all arcs)
     949    ///   ns.reset();
     950    ///   ns.capacityMap(cap).costMap(cost)
     951    ///     .supplyMap(sup).run();
     952    /// \endcode
     953    ///
     954    /// \return <tt>(*this)</tt>
     955    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
     973      return *this;
     974    }
     975
     976    /// @}
     977
     978    /// \name Query Functions
     979    /// The results of the algorithm can be obtained using these
     980    /// functions.\n
     981    /// The \ref run() function must be called before using them.
     982
     983    /// @{
     984
     985    /// \brief Return the total cost of the found flow.
     986    ///
     987    /// This function returns the total cost of the found flow.
     988    /// The complexity of the function is O(e).
     989    ///
     990    /// \note The return type of the function can be specified as a
     991    /// template parameter. For example,
     992    /// \code
     993    ///   ns.totalCost<double>();
     994    /// \endcode
     995    /// It is useful if the total cost cannot be stored in the \c Cost
     996    /// type of the algorithm, which is the default return type of the
     997    /// function.
     998    ///
     999    /// \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];
     1009      }
     1010      return c;
     1011    }
     1012
     1013#ifndef DOXYGEN
     1014    Cost totalCost() const {
     1015      return totalCost<Cost>();
     1016    }
     1017#endif
     1018
     1019    /// \brief Return the flow on the given arc.
     1020    ///
     1021    /// This function returns the flow on the given arc.
     1022    ///
     1023    /// \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.
     1032    ///
     1033    /// \pre \ref run() must be called before using this function.
     1034    const FlowMap& flowMap() const {
     1035      return *_flow_map;
     1036    }
     1037
     1038    /// \brief Return the potential (dual value) of the given node.
     1039    ///
     1040    /// This function returns the potential (dual value) of the
     1041    /// given node.
     1042    ///
     1043    /// \pre \ref run() must be called before using this function.
     1044    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.
     1054    ///
     1055    /// \pre \ref run() must be called before using this function.
     1056    const PotentialMap& potentialMap() const {
     1057      return *_potential_map;
     1058    }
     1059
     1060    /// @}
     1061
     1062  private:
     1063
     1064    // Initialize internal data structures
     1065    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
    6781077      _node_num = countNodes(_graph);
    6791078      _arc_num = countArcs(_graph);
    6801079      int all_node_num = _node_num + 1;
    6811080      int all_arc_num = _arc_num + _node_num;
    682 
     1081      if (_node_num == 0) return false;
     1082
     1083      _arc_ref.resize(_arc_num);
    6831084      _source.resize(all_arc_num);
    6841085      _target.resize(all_arc_num);
    6851086
    686       _lower.resize(all_arc_num);
    687       _upper.resize(all_arc_num);
    6881087      _cap.resize(all_arc_num);
    6891088      _cost.resize(all_arc_num);
    … …  
    7011100      _state.resize(all_arc_num);
    7021101
    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;
    728     }
    729 
    730     /// \name Parameters
    731     /// The parameters of the algorithm can be specified using these
    732     /// functions.
    733 
    734     /// @{
    735 
    736     /// \brief Set the lower bounds on the arcs.
    737     ///
    738     /// This function sets the lower bounds on the arcs.
    739     /// If it is not used before calling \ref run(), the lower bounds
    740     /// will be set to zero on all arcs.
    741     ///
    742     /// \param map An arc map storing the lower bounds.
    743     /// Its \c Value type must be convertible to the \c Value type
    744     /// of the algorithm.
    745     ///
    746     /// \return <tt>(*this)</tt>
    747     template <typename LowerMap>
    748     NetworkSimplex& lowerMap(const LowerMap& map) {
    749       _have_lower = true;
    750       for (ArcIt a(_graph); a != INVALID; ++a) {
    751         _lower[_arc_id[a]] = map[a];
    752       }
    753       return *this;
    754     }
    755 
    756     /// \brief Set the upper bounds (capacities) on the arcs.
    757     ///
    758     /// This function sets the upper bounds (capacities) on the 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).
    762     ///
    763     /// \param map An arc map storing the upper bounds.
    764     /// Its \c Value type must be convertible to the \c Value type
    765     /// of the algorithm.
    766     ///
    767     /// \return <tt>(*this)</tt>
    768     template<typename UpperMap>
    769     NetworkSimplex& upperMap(const UpperMap& map) {
    770       for (ArcIt a(_graph); a != INVALID; ++a) {
    771         _upper[_arc_id[a]] = map[a];
    772       }
    773       return *this;
    774     }
    775 
    776     /// \brief Set the costs of the arcs.
    777     ///
    778     /// This function sets the costs of the arcs.
    779     /// If it is not used before calling \ref run(), the costs
    780     /// will be set to \c 1 on all arcs.
    781     ///
    782     /// \param map An arc map storing the costs.
    783     /// Its \c Value type must be convertible to the \c Cost type
    784     /// of the algorithm.
    785     ///
    786     /// \return <tt>(*this)</tt>
    787     template<typename CostMap>
    788     NetworkSimplex& costMap(const CostMap& map) {
    789       for (ArcIt a(_graph); a != INVALID; ++a) {
    790         _cost[_arc_id[a]] = map[a];
    791       }
    792       return *this;
    793     }
    794 
    795     /// \brief Set the supply values of the nodes.
    796     ///
    797     /// This function sets the supply values of the nodes.
    798     /// If neither this function nor \ref stSupply() is used before
    799     /// calling \ref run(), the supply of each node will be set to zero.
    800     /// (It makes sense only if non-zero lower bounds are given.)
    801     ///
    802     /// \param map A node map storing the supply values.
    803     /// Its \c Value type must be convertible to the \c Value type
    804     /// of the algorithm.
    805     ///
    806     /// \return <tt>(*this)</tt>
    807     template<typename SupplyMap>
    808     NetworkSimplex& supplyMap(const SupplyMap& map) {
    809       for (NodeIt n(_graph); n != INVALID; ++n) {
    810         _supply[_node_id[n]] = map[n];
    811       }
    812       return *this;
    813     }
    814 
    815     /// \brief Set single source and target nodes and a supply value.
    816     ///
    817     /// This function sets a single source node and a single target node
    818     /// and the required flow value.
    819     /// If neither this function nor \ref supplyMap() is used before
    820     /// calling \ref run(), the supply of each node will be set to zero.
    821     /// (It makes sense only if non-zero lower bounds are given.)
    822     ///
    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     ///
    827     /// \param s The source node.
    828     /// \param t The target node.
    829     /// \param k The required amount of flow from node \c s to node \c t
    830     /// (i.e. the supply of \c s and the demand of \c t).
    831     ///
    832     /// \return <tt>(*this)</tt>
    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;
    839       return *this;
    840     }
    841    
    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
    846     /// type will be used.
    847     ///
    848     /// For more information see \ref SupplyType.
    849     ///
    850     /// \return <tt>(*this)</tt>
    851     NetworkSimplex& supplyType(SupplyType supply_type) {
    852       _stype = supply_type;
    853       return *this;
    854     }
    855 
    856     /// @}
    857 
    858     /// \name Execution Control
    859     /// The algorithm can be executed using \ref run().
    860 
    861     /// @{
    862 
    863     /// \brief Run the algorithm.
    864     ///
    865     /// This function runs the algorithm.
    866     /// The paramters can be specified using functions \ref lowerMap(),
    867     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
    868     /// \ref supplyType().
    869     /// For example,
    870     /// \code
    871     ///   NetworkSimplex<ListDigraph> ns(graph);
    872     ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
    873     ///     .supplyMap(sup).run();
    874     /// \endcode
    875     ///
    876     /// This function can be called more than once. All the parameters
    877     /// that have been given are kept for the next call, unless
    878     /// \ref reset() is called, thus only the modified parameters
    879     /// 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.
    882     ///
    883     /// \param pivot_rule The pivot rule that will be used during the
    884     /// algorithm. For more information see \ref PivotRule.
    885     ///
    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);
    898     }
    899 
    900     /// \brief Reset all the parameters that have been given before.
    901     ///
    902     /// This function resets all the paramaters that have been given
    903     /// before using functions \ref lowerMap(), \ref upperMap(),
    904     /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
    905     ///
    906     /// It is useful for multiple run() calls. If this function is not
    907     /// used, all the parameters given before are kept for the next
    908     /// \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.
    911     ///
    912     /// For example,
    913     /// \code
    914     ///   NetworkSimplex<ListDigraph> ns(graph);
    915     ///
    916     ///   // First run
    917     ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
    918     ///     .supplyMap(sup).run();
    919     ///
    920     ///   // Run again with modified cost map (reset() is not called,
    921     ///   // so only the cost map have to be set again)
    922     ///   cost[e] += 100;
    923     ///   ns.costMap(cost).run();
    924     ///
    925     ///   // Run again from scratch using reset()
    926     ///   // (the lower bounds will be set to zero on all arcs)
    927     ///   ns.reset();
    928     ///   ns.upperMap(capacity).costMap(cost)
    929     ///     .supplyMap(sup).run();
    930     /// \endcode
    931     ///
    932     /// \return <tt>(*this)</tt>
    933     NetworkSimplex& reset() {
    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;
    944       return *this;
    945     }
    946 
    947     /// @}
    948 
    949     /// \name Query Functions
    950     /// The results of the algorithm can be obtained using these
    951     /// functions.\n
    952     /// The \ref run() function must be called before using them.
    953 
    954     /// @{
    955 
    956     /// \brief Return the total cost of the found flow.
    957     ///
    958     /// This function returns the total cost of the found flow.
    959     /// Its complexity is O(e).
    960     ///
    961     /// \note The return type of the function can be specified as a
    962     /// template parameter. For example,
    963     /// \code
    964     ///   ns.totalCost<double>();
    965     /// \endcode
    966     /// It is useful if the total cost cannot be stored in the \c Cost
    967     /// type of the algorithm, which is the default return type of the
    968     /// function.
    969     ///
    970     /// \pre \ref run() must be called before using this function.
    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]);
    977       }
    978       return c;
    979     }
    980 
    981 #ifndef DOXYGEN
    982     Cost totalCost() const {
    983       return totalCost<Cost>();
    984     }
    985 #endif
    986 
    987     /// \brief Return the flow on the given arc.
    988     ///
    989     /// This function returns the flow on the given arc.
    990     ///
    991     /// \pre \ref run() must be called before using this function.
    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.
    1001     ///
    1002     /// \pre \ref run() must be called before using this function.
    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       }
    1008     }
    1009 
    1010     /// \brief Return the potential (dual value) of the given node.
    1011     ///
    1012     /// This function returns the potential (dual value) of the
    1013     /// given node.
    1014     ///
    1015     /// \pre \ref run() must be called before using this function.
    1016     Cost potential(const Node& n) const {
    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.
    1026     ///
    1027     /// \pre \ref run() must be called before using this function.
    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       }
    1033     }
    1034 
    1035     /// @}
    1036 
    1037   private:
    1038 
    1039     // Initialize internal data structures
    1040     bool init() {
    1041       if (_node_num == 0) return false;
    1042 
    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) {
     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);
     1119      } 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();
     1135
     1136      // Initialize artifical cost
     1137      Cost art_cost;
     1138      if (std::numeric_limits<Cost>::is_exact) {
     1139        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
     1140      } else {
     1141        art_cost = std::numeric_limits<Cost>::min();
    10531142        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;
     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();
    10571169          } else {
    1058             _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
    1059           }
    1060           _supply[_source[i]] -= c;
    1061           _supply[_target[i]] += c;
     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          }
    10621184        }
    10631185      } else {
    1064         for (int i = 0; i != _arc_num; ++i) {
    1065           _cap[i] = _upper[i];
    1066         }
    1067       }
    1068 
    1069       // Initialize artifical cost
    1070       Cost ART_COST;
    1071       if (std::numeric_limits<Cost>::is_exact) {
    1072         ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
    1073       } else {
    1074         ART_COST = std::numeric_limits<Cost>::min();
    1075         for (int i = 0; i != _arc_num; ++i) {
    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      
     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
    10871216      // Set data for the artificial root node
    10881217      _root = _node_num;
    … …  
    10911220      _thread[_root] = 0;
    10921221      _rev_thread[0] = _root;
    1093       _succ_num[_root] = _node_num + 1;
     1222      _succ_num[_root] = all_node_num;
    10941223      _last_succ[_root] = _root - 1;
    1095       _supply[_root] = -_sum_supply;
    1096       _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
     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      }
    10971285
    10981286      // Add artificial arcs and initialize the spanning tree data structure
    10991287      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
    1100         _parent[u] = _root;
    1101         _pred[u] = e;
    11021288        _thread[u] = u + 1;
    11031289        _rev_thread[u + 1] = u;
    11041290        _succ_num[u] = 1;
    11051291        _last_succ[u] = u;
    1106         _cost[e] = ART_COST;
    1107         _cap[e] = INF;
     1292        _parent[u] = _root;
     1293        _pred[u] = e;
     1294        _cost[e] = art_cost;
     1295        _cap[e] = inf_cap;
    11081296        _state[e] = STATE_TREE;
    1109         if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
     1297        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
    11101298          _flow[e] = _supply[u];
    11111299          _forward[u] = true;
    1112           _pi[u] = -ART_COST + _pi[_root];
     1300          _pi[u] = -art_cost + _pi[_root];
    11131301        } else {
    11141302          _flow[e] = -_supply[u];
    11151303          _forward[u] = false;
    1116           _pi[u] = ART_COST + _pi[_root];
     1304          _pi[u] = art_cost + _pi[_root];
    11171305        }
    11181306      }
    … …  
    11491337      delta = _cap[in_arc];
    11501338      int result = 0;
    1151       Value d;
     1339      Flow d;
    11521340      int e;
    11531341
    … …  
    11551343      for (int u = first; u != join; u = _parent[u]) {
    11561344        e = _pred[u];
    1157         d = _forward[u] ?
    1158           _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
     1345        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
    11591346        if (d < delta) {
    11601347          delta = d;
    … …  
    11661353      for (int u = second; u != join; u = _parent[u]) {
    11671354        e = _pred[u];
    1168         d = _forward[u] ?
    1169           (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
     1355        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
    11701356        if (d <= delta) {
    11711357          delta = d;
    … …  
    11891375      // Augment along the cycle
    11901376      if (delta > 0) {
    1191         Value val = _state[in_arc] * delta;
     1377        Flow val = _state[in_arc] * delta;
    11921378        _flow[in_arc] += val;
    11931379        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
    … …  
    13411527
    13421528    // Execute the algorithm
    1343     ProblemType start(PivotRule pivot_rule) {
     1529    bool start(PivotRule pivot_rule) {
    13441530      // Select the pivot rule implementation
    13451531      switch (pivot_rule) {
    … …  
    13551541          return start<AlteringListPivotRule>();
    13561542      }
    1357       return INFEASIBLE; // avoid warning
     1543      return false;
    13581544    }
    13591545
    13601546    template <typename PivotRuleImpl>
    1361     ProblemType start() {
     1547    bool start() {
    13621548      PivotRuleImpl pivot(*this);
    13631549
    … …  
    13661552        findJoinNode();
    13671553        bool change = findLeavingArc();
    1368         if (delta >= INF) return UNBOUNDED;
    13691554        changeFlow(change);
    13701555        if (change) {
    … …  
    13731558        }
    13741559      }
    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) {
     1560
     1561      // Copy flow values to _flow_map
     1562      if (_plower) {
    13951563        for (int i = 0; i != _arc_num; ++i) {
    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;
     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;
    14061578    }
    14071579

  • lemon/preflow.h

    r641 r611  
    4747
    4848    /// \brief The type of the flow values.
    49     typedef typename CapacityMap::Value Value;
     49    typedef typename CapacityMap::Value Flow;
    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<Value> FlowMap;
     55    typedef typename Digraph::template ArcMap<Flow> FlowMap;
    5656
    5757    /// \brief Instantiates a FlowMap.
    … …  
    8585    ///
    8686    /// The tolerance used by the algorithm to handle inexact computation.
    87     typedef lemon::Tolerance<Value> Tolerance;
     87    typedef lemon::Tolerance<Flow> Tolerance;
    8888
    8989  };
    … …  
    126126    typedef typename Traits::CapacityMap CapacityMap;
    127127    ///The type of the flow values.
    128     typedef typename Traits::Value Value;
     128    typedef typename Traits::Flow Flow;
    129129
    130130    ///The type of the flow map.
    … …  
    152152    bool _local_level;
    153153
    154     typedef typename Digraph::template NodeMap<Value> ExcessMap;
     154    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
    155155    ExcessMap* _excess;
    156156
    … …  
    471471
    472472      for (NodeIt n(_graph); n != INVALID; ++n) {
    473         Value excess = 0;
     473        Flow 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         Value rem = (*_capacity)[e] - (*_flow)[e];
     522        Flow 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         Value rem = (*_flow)[e];
     534        Flow rem = (*_flow)[e];
    535535        if (_tolerance.positive(rem)) {
    536536          Node v = _graph.source(e);
    … …  
    565565
    566566        while (num > 0 && n != INVALID) {
    567           Value excess = (*_excess)[n];
     567          Flow excess = (*_excess)[n];
    568568          int new_level = _level->maxLevel();
    569569
    570570          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
    571             Value rem = (*_capacity)[e] - (*_flow)[e];
     571            Flow 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             Value rem = (*_flow)[e];
     594            Flow 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           Value excess = (*_excess)[n];
     640          Flow excess = (*_excess)[n];
    641641          int new_level = _level->maxLevel();
    642642
    643643          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
    644             Value rem = (*_capacity)[e] - (*_flow)[e];
     644            Flow 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             Value rem = (*_flow)[e];
     667            Flow 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         Value excess = (*_excess)[n];
     781        Flow 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           Value rem = (*_capacity)[e] - (*_flow)[e];
     786          Flow 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           Value rem = (*_flow)[e];
     809          Flow 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     Value flowValue() const {
     900    Flow flowValue() const {
    901901      return (*_excess)[_target];
    902902    }
    903903
    904     /// \brief Returns the flow value on the given arc.
    905     ///
    906     /// Returns the flow value on the given arc. This method can
     904    /// \brief Returns the flow on the given arc.
     905    ///
     906    /// Returns the flow 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     Value flow(const Arc& arc) const {
     911    Flow flow(const Arc& arc) const {
    912912      return (*_flow)[arc];
    913913    }

  • test/min_cost_flow_test.cc

    r642 r615  
    1919#include <iostream>
    2020#include <fstream>
    21 #include <limits>
    2221
    2322#include <lemon/list_graph.h>
    … …  
    3534char test_lgf[] =
    3635  "@nodes\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"               
     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"
    5150  "@arcs\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"
     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"
    7473  "\n"
    7574  "@attributes\n"
    … …  
    7877
    7978
    80 enum SupplyType {
     79enum ProblemType {
    8180  EQ,
    8281  GEQ,
    … …  
    8584
    8685// Check the interface of an MCF algorithm
    87 template <typename GR, typename Value, typename Cost>
     86template <typename GR, typename Flow, typename Cost>
    8887class McfClassConcept
    8988{
    … …  
    9695
    9796      MCF mcf(g);
    98       const MCF& const_mcf = mcf;
    9997
    10098      b = mcf.reset()
    10199             .lowerMap(lower)
    102100             .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)
    106108             .run();
    107 
    108       c = const_mcf.totalCost();
    109       x = const_mcf.template totalCost<double>();
     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>();
    110117      v = const_mcf.flow(a);
    111       c = const_mcf.potential(n);
    112       const_mcf.flowMap(fm);
    113       const_mcf.potentialMap(pm);
     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);
    114123    }
    115124
    116125    typedef typename GR::Node Node;
    117126    typedef typename GR::Arc Arc;
    118     typedef concepts::ReadMap<Node, Value> NM;
    119     typedef concepts::ReadMap<Arc, Value> VAM;
     127    typedef concepts::ReadMap<Node, Flow> NM;
     128    typedef concepts::ReadMap<Arc, Flow> FAM;
    120129    typedef concepts::ReadMap<Arc, Cost> CAM;
    121     typedef concepts::WriteMap<Arc, Value> FlowMap;
    122     typedef concepts::WriteMap<Node, Cost> PotMap;
    123130
    124131    const GR &g;
    125     const VAM &lower;
    126     const VAM &upper;
     132    const FAM &lower;
     133    const FAM &upper;
    127134    const CAM &cost;
    128135    const NM &sup;
    129136    const Node &n;
    130137    const Arc &a;
    131     const Value &k;
    132     FlowMap fm;
    133     PotMap pm;
     138    const Flow &k;
     139    Flow v;
    134140    bool b;
    135     double x;
    136     typename MCF::Value v;
    137     typename MCF::Cost c;
     141
     142    typename MCF::FlowMap &flow;
     143    typename MCF::PotentialMap &pot;
    138144  };
    139145
    … …  
    146152bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
    147153                const SM& supply, const FM& flow,
    148                 SupplyType type = EQ )
     154                ProblemType type = EQ )
    149155{
    150156  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    … …  
    203209template < typename MCF, typename GR,
    204210           typename LM, typename UM,
    205            typename CM, typename SM,
    206            typename PT >
    207 void checkMcf( const MCF& mcf, PT mcf_result,
     211           typename CM, typename SM >
     212void checkMcf( const MCF& mcf, bool mcf_result,
    208213               const GR& gr, const LM& lower, const UM& upper,
    209214               const CM& cost, const SM& supply,
    210                PT result, bool optimal, typename CM::Value total,
     215               bool result, typename CM::Value total,
    211216               const std::string &test_id = "",
    212                SupplyType type = EQ )
     217               ProblemType type = EQ )
    213218{
    214219  check(mcf_result == result, "Wrong result " + test_id);
    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),
     220  if (result) {
     221    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
    221222          "The flow is not feasible " + test_id);
    222223    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
    223     check(checkPotential(gr, lower, upper, cost, supply, flow, pi),
     224    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
     225                         mcf.potentialMap()),
    224226          "Wrong potentials " + test_id);
    225227  }
    … …  
    230232  // Check the interfaces
    231233  {
     234    typedef int Flow;
     235    typedef int Cost;
    232236    typedef concepts::Digraph GR;
    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> >();
     237    checkConcept< McfClassConcept<GR, Flow, Cost>,
     238                  NetworkSimplex<GR, Flow, Cost> >();
    239239  }
    240240
    … …  
    245245  // Read the test digraph
    246246  Digraph 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);
     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);
    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)
    259258    .nodeMap("sup1", s1)
    260259    .nodeMap("sup2", s2)
    … …  
    262261    .nodeMap("sup4", s4)
    263262    .nodeMap("sup5", s5)
    264     .nodeMap("sup6", s6)
    265263    .node("source", v)
    266264    .node("target", w)
    267265    .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;
    308266
    309267  // A. Test NetworkSimplex with the default pivot rule
    … …  
    314272    mcf.upperMap(u).costMap(c);
    315273    checkMcf(mcf, mcf.supplyMap(s1).run(),
    316              gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
     274             gr, l1, u, c, s1, true,  5240, "#A1");
    317275    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
    318              gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
     276             gr, l1, u, c, s2, true,  7620, "#A2");
    319277    mcf.lowerMap(l2);
    320278    checkMcf(mcf, mcf.supplyMap(s1).run(),
    321              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
     279             gr, l2, u, c, s1, true,  5970, "#A3");
    322280    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
    323              gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
     281             gr, l2, u, c, s2, true,  8010, "#A4");
    324282    mcf.reset();
    325283    checkMcf(mcf, mcf.supplyMap(s1).run(),
    326              gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
     284             gr, l1, cu, cc, s1, true,  74, "#A5");
    327285    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
    328              gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
     286             gr, l2, cu, cc, s2, true,  94, "#A6");
    329287    mcf.reset();
    330288    checkMcf(mcf, mcf.run(),
    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");
     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");
    337292
    338293    // Check the GEQ form
    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);
     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);
    343298    checkMcf(mcf, mcf.lowerMap(l2).run(),
    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);
     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);
    348303
    349304    // Check the LEQ form
    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);
     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);
    354309    checkMcf(mcf, mcf.lowerMap(l2).run(),
    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");
     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);
    370314  }
    371315
    … …  
    373317  {
    374318    NetworkSimplex<Digraph> mcf(gr);
    375     mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
     319    mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
    376320
    377321    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
    378              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
     322             gr, l2, u, c, s1, true,  5970, "#B1");
    379323    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
    380              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
     324             gr, l2, u, c, s1, true,  5970, "#B2");
    381325    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
    382              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
     326             gr, l2, u, c, s1, true,  5970, "#B3");
    383327    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
    384              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
     328             gr, l2, u, c, s1, true,  5970, "#B4");
    385329    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
    386              gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
     330             gr, l2, u, c, s1, true,  5970, "#B5");
    387331  }
    388332

  • tools/dimacs-solver.cc

    r644 r627  
    120120  ti.restart();
    121121  NetworkSimplex<Digraph, Value> ns(g);
    122   ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup);
    123   if (sum_sup > 0) ns.supplyType(ns.LEQ);
     122  ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
     123  if (sum_sup > 0) ns.problemType(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.