COIN-OR::LEMON - Graph Library

About

Downloads

Documentation

Development

Issue Tracker

Browse Tickets

Custom Query

Roadmap

Browse Source

Local Info

Mailing Lists

Search

Ticket #270: ns-negative-6c408d864fa1.patch

File ns-negative-6c408d864fa1.patch, 43.8 KB (added by Peter Kovacs, 15 years ago)

  • doc/groups.dox 

    # HG changeset patch
# User Peter Kovacs <kpeter@inf.elte.hu>
# Date 1240967724 -7200
# Node ID 6c408d864fa1066df1ba3f323b9396d940c9cc19
# Parent  58357e986a08f59fd3fe0db28468b69517950a31
Support negative costs and bounds in NetworkSimplex (#270)

  * The interface is reworked to support negative costs and bounds.
    - ProblemType and problemType() are renamed to
      SupplyType and supplyType(), see also #234.
    - ProblemType type is introduced similarly to the LP interface.
    - 'bool run()' is replaced by 'ProblemType run()' to handle
      unbounded problem instances, as well.
    - Add INF public member constant similarly to the LP interface.
  * Remove capacityMap() and boundMaps(), see also #266.
  * Update the problem definition in the MCF module.
  * Remove the usage of Circulation (and adaptors) for checking feasibility.
    Check feasibility by examining the artifical arcs instead (after solving
    the problem).
  * Additional check for unbounded negative cycles found during the
    algorithm (it is possible now, since negative costs are allowed).
  * Fix in the constructor (the value types needn't be integer any more),
    see also #254.
  * Improve and extend the doc.
  * Rework the test file and add test cases for negative costs and bounds.

diff --git a/doc/groups.dox b/doc/groups.dox
    a b  
    352352minimum total cost from a set of supply nodes to a set of demand nodes
    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
    368368\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    … …  
    404404
    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
    410410conditions hold.
    … …  
    413413   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
    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
    427427LEMON contains several algorithms for solving minimum cost flow problems.

  • lemon/network_simplex.h 

    diff --git a/lemon/network_simplex.h b/lemon/network_simplex.h
    a b  
    3030
    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 {
    3835
    … …  
    5047  ///
    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.
    5759  /// \tparam F The value type used for flow amounts, capacity bounds
    … …  
    8890
    8991  public:
    9092
    91     /// \brief Enum type for selecting the pivot rule.
     93    /// \brief Problem type constants for the \c run() function.
    9294    ///
    93     /// Enum type for selecting the pivot rule for the \ref run()
     95    /// Enum type containing the problem type constants that can be
     96    /// returned by the \ref run() function of the algorithm.
     97    enum ProblemType {
     98      /// The problem has no feasible solution (flow).
     99      INFEASIBLE,
     100      /// The problem has optimal solution (i.e. it is feasible and
     101      /// bounded), and the algorithm has found optimal flow and node
     102      /// potentials (primal and dual solutions).
     103      OPTIMAL,
     104      /// The objective function of the problem is unbounded, i.e.
     105      /// there is a directed cycle having negative total cost and
     106      /// infinite upper bound.
     107      UNBOUNDED
     108    };
     109   
     110    /// \brief Constants for selecting the type of the supply constraints.
     111    ///
     112    /// Enum type containing constants for selecting the supply type,
     113    /// i.e. the direction of the inequalities in the supply/demand
     114    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
     115    ///
     116    /// The default supply type is \c GEQ, since this form is supported
     117    /// by other minimum cost flow algorithms and the \ref Circulation
     118    /// algorithm, as well.
     119    /// The \c LEQ problem type can be selected using the \ref supplyType()
    94120    /// function.
    95121    ///
     122    /// Note that the equality form is a special case of both supply types.
     123    enum SupplyType {
     124
     125      /// This option means that there are <em>"greater or equal"</em>
     126      /// supply/demand constraints in the definition, i.e. the exact
     127      /// formulation of the problem is the following.
     128      /**
     129          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
     130          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
     131              sup(u) \quad \forall u\in V \f]
     132          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
     133      */
     134      /// It means that the total demand must be greater or equal to the
     135      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
     136      /// negative) and all the supplies have to be carried out from
     137      /// the supply nodes, but there could be demands that are not
     138      /// satisfied.
     139      GEQ,
     140      /// It is just an alias for the \c GEQ option.
     141      CARRY_SUPPLIES = GEQ,
     142
     143      /// This option means that there are <em>"less or equal"</em>
     144      /// supply/demand constraints in the definition, i.e. the exact
     145      /// formulation of the problem is the following.
     146      /**
     147          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
     148          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
     149              sup(u) \quad \forall u\in V \f]
     150          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
     151      */
     152      /// It means that the total demand must be less or equal to the
     153      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
     154      /// positive) and all the demands have to be satisfied, but there
     155      /// could be supplies that are not carried out from the supply
     156      /// nodes.
     157      LEQ,
     158      /// It is just an alias for the \c LEQ option.
     159      SATISFY_DEMANDS = LEQ
     160    };
     161   
     162    /// \brief Constants for selecting the pivot rule.
     163    ///
     164    /// Enum type containing constants for selecting the pivot rule for
     165    /// the \ref run() function.
     166    ///
    96167    /// \ref NetworkSimplex provides five different pivot rule
    97168    /// implementations that significantly affect the running time
    98169    /// of the algorithm.
    … …  
    131202      ALTERING_LIST
    132203    };
    133204   
    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
    141     /// 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()
    144     /// function.
    145     ///
    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.
    152       /**
    153           \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    154           \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
    155               sup(u) \quad \forall u\in V \f]
    156           \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    157       */
    158       /// It means that the total demand must be greater or equal to the
    159       /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
    160       /// negative) and all the supplies have to be carried out from
    161       /// the supply nodes, but there could be demands that are not
    162       /// satisfied.
    163       GEQ,
    164       /// It is just an alias for the \c GEQ option.
    165       CARRY_SUPPLIES = GEQ,
    166 
    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.
    170       /**
    171           \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    172           \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
    173               sup(u) \quad \forall u\in V \f]
    174           \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    175       */
    176       /// It means that the total demand must be less or equal to the
    177       /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
    178       /// positive) and all the demands have to be satisfied, but there
    179       /// could be supplies that are not carried out from the supply
    180       /// nodes.
    181       LEQ,
    182       /// It is just an alias for the \c LEQ option.
    183       SATISFY_DEMANDS = LEQ
    184     };
    185 
    186205  private:
    187206
    188207    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    … …  
    220239    bool _pstsup;
    221240    Node _psource, _ptarget;
    222241    Flow _pstflow;
    223     ProblemType _ptype;
     242    SupplyType _stype;
     243   
     244    Flow _sum_supply;
    224245
    225246    // Result maps
    226247    FlowMap *_flow_map;
    … …  
    259280    int stem, par_stem, new_stem;
    260281    Flow delta;
    261282
     283  public:
     284 
     285    /// \brief Constant for infinite upper bounds (capacities).
     286    ///
     287    /// Constant for infinite upper bounds (capacities).
     288    /// It is \c std::numeric_limits<Flow>::infinity() if available,
     289    /// \c std::numeric_limits<Flow>::max() otherwise.
     290    const Flow INF;
     291
    262292  private:
    263293
    264294    // Implementation of the First Eligible pivot rule
    … …  
    661691    NetworkSimplex(const GR& graph) :
    662692      _graph(graph),
    663693      _plower(NULL), _pupper(NULL), _pcost(NULL),
    664       _psupply(NULL), _pstsup(false), _ptype(GEQ),
     694      _psupply(NULL), _pstsup(false), _stype(GEQ),
    665695      _flow_map(NULL), _potential_map(NULL),
    666696      _local_flow(false), _local_potential(false),
    667       _node_id(graph)
     697      _node_id(graph),
     698      INF(std::numeric_limits<Flow>::has_infinity ?
     699          std::numeric_limits<Flow>::infinity() :
     700          std::numeric_limits<Flow>::max())
    668701    {
    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");
     702      // Check the value types
     703      LEMON_ASSERT(std::numeric_limits<Flow>::is_signed,
     704        "The flow type of NetworkSimplex must be signed");
     705      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
     706        "The cost type of NetworkSimplex must be signed");
    675707    }
    676708
    677709    /// Destructor.
    … …  
    689721    /// \brief Set the lower bounds on the arcs.
    690722    ///
    691723    /// 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.
     724    /// If it is not used before calling \ref run(), the lower bounds
     725    /// will be set to zero on all arcs.
    695726    ///
    696727    /// \param map An arc map storing the lower bounds.
    697728    /// Its \c Value type must be convertible to the \c Flow type
    698729    /// of the algorithm.
    699730    ///
    700731    /// \return <tt>(*this)</tt>
    701     template <typename LOWER>
    702     NetworkSimplex& lowerMap(const LOWER& map) {
     732    template <typename LowerMap>
     733    NetworkSimplex& lowerMap(const LowerMap& map) {
    703734      delete _plower;
    704735      _plower = new FlowArcMap(_graph);
    705736      for (ArcIt a(_graph); a != INVALID; ++a) {
    … …  
    711742    /// \brief Set the upper bounds (capacities) on the arcs.
    712743    ///
    713744    /// 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.
     745    /// If it is not used before calling \ref run(), the upper bounds
     746    /// will be set to \ref INF on all arcs (i.e. the flow value will be
     747    /// unbounded from above on each arc).
    718748    ///
    719749    /// \param map An arc map storing the upper bounds.
    720750    /// Its \c Value type must be convertible to the \c Flow type
    721751    /// of the algorithm.
    722752    ///
    723753    /// \return <tt>(*this)</tt>
    724     template<typename UPPER>
    725     NetworkSimplex& upperMap(const UPPER& map) {
     754    template<typename UpperMap>
     755    NetworkSimplex& upperMap(const UpperMap& map) {
    726756      delete _pupper;
    727757      _pupper = new FlowArcMap(_graph);
    728758      for (ArcIt a(_graph); a != INVALID; ++a) {
    … …  
    731761      return *this;
    732762    }
    733763
    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 
    771764    /// \brief Set the costs of the arcs.
    772765    ///
    773766    /// This function sets the costs of the arcs.
    … …  
    779772    /// of the algorithm.
    780773    ///
    781774    /// \return <tt>(*this)</tt>
    782     template<typename COST>
    783     NetworkSimplex& costMap(const COST& map) {
     775    template<typename CostMap>
     776    NetworkSimplex& costMap(const CostMap& map) {
    784777      delete _pcost;
    785778      _pcost = new CostArcMap(_graph);
    786779      for (ArcIt a(_graph); a != INVALID; ++a) {
    … …  
    801794    /// of the algorithm.
    802795    ///
    803796    /// \return <tt>(*this)</tt>
    804     template<typename SUP>
    805     NetworkSimplex& supplyMap(const SUP& map) {
     797    template<typename SupplyMap>
     798    NetworkSimplex& supplyMap(const SupplyMap& map) {
    806799      delete _psupply;
    807800      _pstsup = false;
    808801      _psupply = new FlowNodeMap(_graph);
    … …  
    820813    /// calling \ref run(), the supply of each node will be set to zero.
    821814    /// (It makes sense only if non-zero lower bounds are given.)
    822815    ///
     816    /// Using this function has the same effect as using \ref supplyMap()
     817    /// with such a map in which \c k is assigned to \c s, \c -k is
     818    /// assigned to \c t and all other nodes have zero supply value.
     819    ///
    823820    /// \param s The source node.
    824821    /// \param t The target node.
    825822    /// \param k The required amount of flow from node \c s to node \c t
    … …  
    836833      return *this;
    837834    }
    838835   
    839     /// \brief Set the problem type.
     836    /// \brief Set the type of the supply constraints.
    840837    ///
    841     /// This function sets the problem type for the algorithm.
    842     /// If it is not used before calling \ref run(), the \ref GEQ problem
     838    /// This function sets the type of the supply/demand constraints.
     839    /// If it is not used before calling \ref run(), the \ref GEQ supply
    843840    /// type will be used.
    844841    ///
    845     /// For more information see \ref ProblemType.
     842    /// For more information see \ref SupplyType.
    846843    ///
    847844    /// \return <tt>(*this)</tt>
    848     NetworkSimplex& problemType(ProblemType problem_type) {
    849       _ptype = problem_type;
     845    NetworkSimplex& supplyType(SupplyType supply_type) {
     846      _stype = supply_type;
    850847      return *this;
    851848    }
    852849
    … …  
    896893    ///
    897894    /// This function runs the algorithm.
    898895    /// 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().
     896    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
     897    /// \ref supplyType(), \ref flowMap() and \ref potentialMap().
    902898    /// For example,
    903899    /// \code
    904900    ///   NetworkSimplex<ListDigraph> ns(graph);
    905     ///   ns.boundMaps(lower, upper).costMap(cost)
     901    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
    906902    ///     .supplyMap(sup).run();
    907903    /// \endcode
    908904    ///
    … …  
    914910    /// \param pivot_rule The pivot rule that will be used during the
    915911    /// algorithm. For more information see \ref PivotRule.
    916912    ///
    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);
     913    /// \return \c INFEASIBLE if no feasible flow exists,
     914    /// \n \c OPTIMAL if the problem has optimal solution
     915    /// (i.e. it is feasible and bounded), and the algorithm has found
     916    /// optimal flow and node potentials (primal and dual solutions),
     917    /// \n \c UNBOUNDED if the objective function of the problem is
     918    /// unbounded, i.e. there is a directed cycle having negative total
     919    /// cost and infinite upper bound.
     920    ///
     921    /// \see ProblemType, PivotRule
     922    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
     923      if (!init()) return INFEASIBLE;
     924      return start(pivot_rule);
    920925    }
    921926
    922927    /// \brief Reset all the parameters that have been given before.
    923928    ///
    924929    /// This function resets all the paramaters that have been given
    925930    /// before using functions \ref lowerMap(), \ref upperMap(),
    926     /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
    927     /// \ref supplyMap(), \ref stSupply(), \ref problemType(),
     931    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(),
    928932    /// \ref flowMap() and \ref potentialMap().
    929933    ///
    930934    /// It is useful for multiple run() calls. If this function is not
    … …  
    936940    ///   NetworkSimplex<ListDigraph> ns(graph);
    937941    ///
    938942    ///   // First run
    939     ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
     943    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
    940944    ///     .supplyMap(sup).run();
    941945    ///
    942946    ///   // Run again with modified cost map (reset() is not called,
    … …  
    947951    ///   // Run again from scratch using reset()
    948952    ///   // (the lower bounds will be set to zero on all arcs)
    949953    ///   ns.reset();
    950     ///   ns.capacityMap(cap).costMap(cost)
     954    ///   ns.upperMap(capacity).costMap(cost)
    951955    ///     .supplyMap(sup).run();
    952956    /// \endcode
    953957    ///
    … …  
    962966      _pcost = NULL;
    963967      _psupply = NULL;
    964968      _pstsup = false;
    965       _ptype = GEQ;
     969      _stype = GEQ;
    966970      if (_local_flow) delete _flow_map;
    967971      if (_local_potential) delete _potential_map;
    968972      _flow_map = NULL;
    … …  
    985989    /// \brief Return the total cost of the found flow.
    986990    ///
    987991    /// This function returns the total cost of the found flow.
    988     /// The complexity of the function is O(e).
     992    /// Its complexity is O(e).
    989993    ///
    990994    /// \note The return type of the function can be specified as a
    991995    /// template parameter. For example,
    … …  
    9971001    /// function.
    9981002    ///
    9991003    /// \pre \ref run() must be called before using this function.
    1000     template <typename Num>
    1001     Num totalCost() const {
    1002       Num c = 0;
     1004    template <typename Value>
     1005    Value totalCost() const {
     1006      Value c = 0;
    10031007      if (_pcost) {
    10041008        for (ArcIt e(_graph); e != INVALID; ++e)
    10051009          c += (*_flow_map)[e] * (*_pcost)[e];
    … …  
    10501054    ///
    10511055    /// This function returns a const reference to a node map storing
    10521056    /// the found potentials, which form the dual solution of the
    1053     /// \ref min_cost_flow "minimum cost flow" problem.
     1057    /// \ref min_cost_flow "minimum cost flow problem".
    10541058    ///
    10551059    /// \pre \ref run() must be called before using this function.
    10561060    const PotentialMap& potentialMap() const {
    … …  
    11011105
    11021106      // Initialize node related data
    11031107      bool valid_supply = true;
    1104       Flow sum_supply = 0;
     1108      _sum_supply = 0;
    11051109      if (!_pstsup && !_psupply) {
    11061110        _pstsup = true;
    11071111        _psource = _ptarget = NodeIt(_graph);
    … …  
    11121116        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
    11131117          _node_id[n] = i;
    11141118          _supply[i] = (*_psupply)[n];
    1115           sum_supply += _supply[i];
     1119          _sum_supply += _supply[i];
    11161120        }
    1117         valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
    1118                        (_ptype == LEQ && sum_supply >= 0);
     1121        valid_supply = (_stype == GEQ && _sum_supply <= 0) ||
     1122                       (_stype == LEQ && _sum_supply >= 0);
    11191123      } else {
    11201124        int i = 0;
    11211125        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
    … …  
    11271131      }
    11281132      if (!valid_supply) return false;
    11291133
    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 
    11361134      // Initialize artifical cost
    1137       Cost art_cost;
     1135      Cost ART_COST;
    11381136      if (std::numeric_limits<Cost>::is_exact) {
    1139         art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
     1137        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
    11401138      } else {
    1141         art_cost = std::numeric_limits<Cost>::min();
     1139        ART_COST = std::numeric_limits<Cost>::min();
    11421140        for (int i = 0; i != _arc_num; ++i) {
    1143           if (_cost[i] > art_cost) art_cost = _cost[i];
     1141          if (_cost[i] > ART_COST) ART_COST = _cost[i];
    11441142        }
    1145         art_cost = (art_cost + 1) * _node_num;
     1143        ART_COST = (ART_COST + 1) * _node_num;
    11461144      }
    11471145
    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 
    12161146      // Set data for the artificial root node
    12171147      _root = _node_num;
    12181148      _parent[_root] = -1;
    … …  
    12211151      _rev_thread[0] = _root;
    12221152      _succ_num[_root] = all_node_num;
    12231153      _last_succ[_root] = _root - 1;
    1224       _supply[_root] = -sum_supply;
    1225       if (sum_supply < 0) {
    1226         _pi[_root] = -art_cost;
     1154      _supply[_root] = -_sum_supply;
     1155      if (_sum_supply < 0) {
     1156        _pi[_root] = -ART_COST;
    12271157      } else {
    1228         _pi[_root] = art_cost;
     1158        _pi[_root] = ART_COST;
    12291159      }
    12301160
    12311161      // Store the arcs in a mixed order
    … …  
    12601190            _cap[i] = (*_pupper)[_arc_ref[i]];
    12611191        } else {
    12621192          for (int i = 0; i != _arc_num; ++i)
    1263             _cap[i] = inf_cap;
     1193            _cap[i] = INF;
    12641194        }
    12651195        if (_pcost) {
    12661196          for (int i = 0; i != _arc_num; ++i)
    … …  
    12751205      if (_plower) {
    12761206        for (int i = 0; i != _arc_num; ++i) {
    12771207          Flow c = (*_plower)[_arc_ref[i]];
    1278           if (c != 0) {
    1279             _cap[i] -= c;
     1208          if (c > 0) {
     1209            if (_cap[i] < INF) _cap[i] -= c;
     1210            _supply[_source[i]] -= c;
     1211            _supply[_target[i]] += c;
     1212          }
     1213          else if (c < 0) {
     1214            if (_cap[i] < INF + c) {
     1215              _cap[i] -= c;
     1216            } else {
     1217              _cap[i] = INF;
     1218            }
    12801219            _supply[_source[i]] -= c;
    12811220            _supply[_target[i]] += c;
    12821221          }
    … …  
    12911230        _last_succ[u] = u;
    12921231        _parent[u] = _root;
    12931232        _pred[u] = e;
    1294         _cost[e] = art_cost;
    1295         _cap[e] = inf_cap;
     1233        _cost[e] = ART_COST;
     1234        _cap[e] = INF;
    12961235        _state[e] = STATE_TREE;
    1297         if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
     1236        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
    12981237          _flow[e] = _supply[u];
    12991238          _forward[u] = true;
    1300           _pi[u] = -art_cost + _pi[_root];
     1239          _pi[u] = -ART_COST + _pi[_root];
    13011240        } else {
    13021241          _flow[e] = -_supply[u];
    13031242          _forward[u] = false;
    1304           _pi[u] = art_cost + _pi[_root];
     1243          _pi[u] = ART_COST + _pi[_root];
    13051244        }
    13061245      }
    13071246
    … …  
    13421281      // Search the cycle along the path form the first node to the root
    13431282      for (int u = first; u != join; u = _parent[u]) {
    13441283        e = _pred[u];
    1345         d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
     1284        d = _forward[u] ?
     1285          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
    13461286        if (d < delta) {
    13471287          delta = d;
    13481288          u_out = u;
    … …  
    13521292      // Search the cycle along the path form the second node to the root
    13531293      for (int u = second; u != join; u = _parent[u]) {
    13541294        e = _pred[u];
    1355         d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
     1295        d = _forward[u] ?
     1296          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
    13561297        if (d <= delta) {
    13571298          delta = d;
    13581299          u_out = u;
    … …  
    15261467    }
    15271468
    15281469    // Execute the algorithm
    1529     bool start(PivotRule pivot_rule) {
     1470    ProblemType start(PivotRule pivot_rule) {
    15301471      // Select the pivot rule implementation
    15311472      switch (pivot_rule) {
    15321473        case FIRST_ELIGIBLE:
    … …  
    15401481        case ALTERING_LIST:
    15411482          return start<AlteringListPivotRule>();
    15421483      }
    1543       return false;
     1484      return INFEASIBLE; // avoid warning
    15441485    }
    15451486
    15461487    template <typename PivotRuleImpl>
    1547     bool start() {
     1488    ProblemType start() {
    15481489      PivotRuleImpl pivot(*this);
    15491490
    15501491      // Execute the Network Simplex algorithm
    15511492      while (pivot.findEnteringArc()) {
    15521493        findJoinNode();
    15531494        bool change = findLeavingArc();
     1495        if (delta >= INF) return UNBOUNDED;
    15541496        changeFlow(change);
    15551497        if (change) {
    15561498          updateTreeStructure();
    15571499          updatePotential();
    15581500        }
    15591501      }
     1502     
     1503      // Check feasibility
     1504      if (_sum_supply < 0) {
     1505        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
     1506          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
     1507        }
     1508      }
     1509      else if (_sum_supply > 0) {
     1510        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
     1511          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
     1512        }
     1513      }
     1514      else {
     1515        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
     1516          if (_flow[e] != 0) return INFEASIBLE;
     1517        }
     1518      }
    15601519
    15611520      // Copy flow values to _flow_map
    15621521      if (_plower) {
    … …  
    15741533        _potential_map->set(n, _pi[_node_id[n]]);
    15751534      }
    15761535
    1577       return true;
     1536      return OPTIMAL;
    15781537    }
    15791538
    15801539  }; //class NetworkSimplex

  • test/min_cost_flow_test.cc 

    diff --git a/test/min_cost_flow_test.cc b/test/min_cost_flow_test.cc
    a b  
    1818
    1919#include <iostream>
    2020#include <fstream>
     21#include <limits>
    2122
    2223#include <lemon/list_graph.h>
    2324#include <lemon/lgf_reader.h>
    … …  
    3334
    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"
    7576  "source 1\n"
    7677  "target 12\n";
    7778
    7879
    79 enum ProblemType {
     80enum SupplyType {
    8081  EQ,
    8182  GEQ,
    8283  LEQ
    … …  
    9899      b = mcf.reset()
    99100             .lowerMap(lower)
    100101             .upperMap(upper)
    101              .capacityMap(upper)
    102              .boundMaps(lower, upper)
    103102             .costMap(cost)
    104103             .supplyMap(sup)
    105104             .stSupply(n, n, k)
    … …  
    112111      const typename MCF::FlowMap &fm = const_mcf.flowMap();
    113112      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
    114113
    115       v = const_mcf.totalCost();
     114      c = const_mcf.totalCost();
    116115      double x = const_mcf.template totalCost<double>();
    117116      v = const_mcf.flow(a);
    118       v = const_mcf.potential(n);
     117      c = const_mcf.potential(n);
     118     
     119      v = const_mcf.INF;
    119120
    120121      ignore_unused_variable_warning(fm);
    121122      ignore_unused_variable_warning(pm);
    … …  
    137138    const Arc &a;
    138139    const Flow &k;
    139140    Flow v;
     141    Cost c;
    140142    bool b;
    141143
    142144    typename MCF::FlowMap &flow;
    … …  
    151153           typename SM, typename FM >
    152154bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
    153155                const SM& supply, const FM& flow,
    154                 ProblemType type = EQ )
     156                SupplyType type = EQ )
    155157{
    156158  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    157159
    … …  
    208210// Run a minimum cost flow algorithm and check the results
    209211template < typename MCF, typename GR,
    210212           typename LM, typename UM,
    211            typename CM, typename SM >
    212 void checkMcf( const MCF& mcf, bool mcf_result,
     213           typename CM, typename SM,
     214           typename PT >
     215void checkMcf( const MCF& mcf, PT mcf_result,
    213216               const GR& gr, const LM& lower, const UM& upper,
    214217               const CM& cost, const SM& supply,
    215                bool result, typename CM::Value total,
     218               PT result, bool optimal, typename CM::Value total,
    216219               const std::string &test_id = "",
    217                ProblemType type = EQ )
     220               SupplyType type = EQ )
    218221{
    219222  check(mcf_result == result, "Wrong result " + test_id);
    220   if (result) {
     223  if (optimal) {
    221224    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
    222225          "The flow is not feasible " + test_id);
    223226    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
    … …  
    244247
    245248  // Read the test digraph
    246249  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);
     250  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
     251  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
    249252  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
    250253  Node v, w;
    251254
    … …  
    255258    .arcMap("cap", u)
    256259    .arcMap("low1", l1)
    257260    .arcMap("low2", l2)
     261    .arcMap("low3", l3)
    258262    .nodeMap("sup1", s1)
    259263    .nodeMap("sup2", s2)
    260264    .nodeMap("sup3", s3)
    261265    .nodeMap("sup4", s4)
    262266    .nodeMap("sup5", s5)
     267    .nodeMap("sup6", s6)
    263268    .node("source", v)
    264269    .node("target", w)
    265270    .run();
     271 
     272  // Build a test digraph for testing negative costs
     273  Digraph ngr;
     274  Node n1 = ngr.addNode();
     275  Node n2 = ngr.addNode();
     276  Node n3 = ngr.addNode();
     277  Node n4 = ngr.addNode();
     278  Node n5 = ngr.addNode();
     279  Node n6 = ngr.addNode();
     280  Node n7 = ngr.addNode();
     281 
     282  Arc a1 = ngr.addArc(n1, n2);
     283  Arc a2 = ngr.addArc(n1, n3);
     284  Arc a3 = ngr.addArc(n2, n4);
     285  Arc a4 = ngr.addArc(n3, n4);
     286  Arc a5 = ngr.addArc(n3, n2);
     287  Arc a6 = ngr.addArc(n5, n3);
     288  Arc a7 = ngr.addArc(n5, n6);
     289  Arc a8 = ngr.addArc(n6, n7);
     290  Arc a9 = ngr.addArc(n7, n5);
     291 
     292  Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
     293  ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
     294  Digraph::NodeMap<int> ns(ngr, 0);
     295 
     296  nl2[a7] =  1000;
     297  nl2[a8] = -1000;
     298 
     299  ns[n1] =  100;
     300  ns[n4] = -100;
     301 
     302  nc[a1] =  100;
     303  nc[a2] =   30;
     304  nc[a3] =   20;
     305  nc[a4] =   80;
     306  nc[a5] =   50;
     307  nc[a6] =   10;
     308  nc[a7] =   80;
     309  nc[a8] =   30;
     310  nc[a9] = -120;
    266311
    267312  // A. Test NetworkSimplex with the default pivot rule
    268313  {
    … …  
    271316    // Check the equality form
    272317    mcf.upperMap(u).costMap(c);
    273318    checkMcf(mcf, mcf.supplyMap(s1).run(),
    274              gr, l1, u, c, s1, true,  5240, "#A1");
     319             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
    275320    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
    276              gr, l1, u, c, s2, true,  7620, "#A2");
     321             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
    277322    mcf.lowerMap(l2);
    278323    checkMcf(mcf, mcf.supplyMap(s1).run(),
    279              gr, l2, u, c, s1, true,  5970, "#A3");
     324             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
    280325    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
    281              gr, l2, u, c, s2, true,  8010, "#A4");
     326             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
    282327    mcf.reset();
    283328    checkMcf(mcf, mcf.supplyMap(s1).run(),
    284              gr, l1, cu, cc, s1, true,  74, "#A5");
     329             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
    285330    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
    286              gr, l2, cu, cc, s2, true,  94, "#A6");
     331             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
    287332    mcf.reset();
    288333    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");
     334             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
     335    checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
     336             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
     337    mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
     338    checkMcf(mcf, mcf.run(),
     339             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
    292340
    293341    // Check the GEQ form
    294     mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
     342    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
    295343    checkMcf(mcf, mcf.run(),
    296              gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
    297     mcf.problemType(mcf.GEQ);
     344             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
     345    mcf.supplyType(mcf.GEQ);
    298346    checkMcf(mcf, mcf.lowerMap(l2).run(),
    299              gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
    300     mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
     347             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
     348    mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
    301349    checkMcf(mcf, mcf.run(),
    302              gr, l2, u, c, s5, false,    0, "#A11", GEQ);
     350             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
    303351
    304352    // Check the LEQ form
    305     mcf.reset().problemType(mcf.LEQ);
    306     mcf.upperMap(u).costMap(c).supplyMap(s5);
     353    mcf.reset().supplyType(mcf.LEQ);
     354    mcf.upperMap(u).costMap(c).supplyMap(s6);
    307355    checkMcf(mcf, mcf.run(),
    308              gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
     356             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
    309357    checkMcf(mcf, mcf.lowerMap(l2).run(),
    310              gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
    311     mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
     358             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
     359    mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
    312360    checkMcf(mcf, mcf.run(),
    313              gr, l2, u, c, s4, false,    0, "#A14", LEQ);
     361             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
     362
     363    // Check negative costs
     364    NetworkSimplex<Digraph> nmcf(ngr);
     365    nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
     366    checkMcf(nmcf, nmcf.run(),
     367      ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
     368    checkMcf(nmcf, nmcf.upperMap(nu2).run(),
     369      ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
     370    nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
     371    checkMcf(nmcf, nmcf.run(),
     372      ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
    314373  }
    315374
    316375  // B. Test NetworkSimplex with each pivot rule
    317376  {
    318377    NetworkSimplex<Digraph> mcf(gr);
    319     mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
     378    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
    320379
    321380    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
    322              gr, l2, u, c, s1, true,  5970, "#B1");
     381             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
    323382    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
    324              gr, l2, u, c, s1, true,  5970, "#B2");
     383             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
    325384    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
    326              gr, l2, u, c, s1, true,  5970, "#B3");
     385             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
    327386    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
    328              gr, l2, u, c, s1, true,  5970, "#B4");
     387             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
    329388    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
    330              gr, l2, u, c, s1, true,  5970, "#B5");
     389             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
    331390  }
    332391
    333392  return 0;

  • tools/dimacs-solver.cc 

    diff --git a/tools/dimacs-solver.cc b/tools/dimacs-solver.cc
    a b  
    119119
    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();
    126126  bool res = ns.run();