COIN-OR::LEMON - Graph Library

About

Downloads

Documentation

Development

Issue Tracker

Browse Tickets

Custom Query

Roadmap

Browse Source

Local Info

Mailing Lists

Search

Changeset 687:6c408d864fa1 in lemon for lemon/network_simplex.h


Ignore:
Timestamp:
04/29/09 03:15:24 (16 years ago)
Author:
Peter Kovacs <kpeter@…>
Branch:
default
Phase:
public
Message:

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.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • lemon/network_simplex.h

    r665 r687  
    3131#include <lemon/core.h>
    3232#include <lemon/math.h>
    33 #include <lemon/maps.h>
    34 #include <lemon/circulation.h>
    35 #include <lemon/adaptors.h>
    3633
    3734namespace lemon {
    … …  
    5148  /// In general this class is the fastest implementation available
    5249  /// in LEMON for the minimum cost flow problem.
    53   /// Moreover it supports both direction of the supply/demand inequality
    54   /// constraints. For more information see \ref ProblemType.
     50  /// Moreover it supports both directions of the supply/demand inequality
     51  /// constraints. For more information see \ref SupplyType.
     52  ///
     53  /// Most of the parameters of the problem (except for the digraph)
     54  /// can be given using separate functions, and the algorithm can be
     55  /// executed using the \ref run() function. If some parameters are not
     56  /// specified, then default values will be used.
    5557  ///
    5658  /// \tparam GR The digraph type the algorithm runs on.
    … …  
    8991  public:
    9092
    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
     93    /// \brief Problem type constants for the \c run() function.
     94    ///
     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
    132108    };
    133109   
    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
     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
    141117    /// 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()
     118    /// algorithm, as well.
     119    /// The \c LEQ problem type can be selected using the \ref supplyType()
    144120    /// function.
    145121    ///
    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.
     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.
    152128      /**
    153129          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    … …  
    165141      CARRY_SUPPLIES = GEQ,
    166142
    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.
     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.
    170146      /**
    171147          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    … …  
    183159      SATISFY_DEMANDS = LEQ
    184160    };
    185 
     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    ///
     167    /// \ref NetworkSimplex provides five different pivot rule
     168    /// implementations that significantly affect the running time
     169    /// of the algorithm.
     170    /// By default \ref BLOCK_SEARCH "Block Search" is used, which
     171    /// proved to be the most efficient and the most robust on various
     172    /// test inputs according to our benchmark tests.
     173    /// However another pivot rule can be selected using the \ref run()
     174    /// function with the proper parameter.
     175    enum PivotRule {
     176
     177      /// The First Eligible pivot rule.
     178      /// The next eligible arc is selected in a wraparound fashion
     179      /// in every iteration.
     180      FIRST_ELIGIBLE,
     181
     182      /// The Best Eligible pivot rule.
     183      /// The best eligible arc is selected in every iteration.
     184      BEST_ELIGIBLE,
     185
     186      /// The Block Search pivot rule.
     187      /// A specified number of arcs are examined in every iteration
     188      /// in a wraparound fashion and the best eligible arc is selected
     189      /// from this block.
     190      BLOCK_SEARCH,
     191
     192      /// The Candidate List pivot rule.
     193      /// In a major iteration a candidate list is built from eligible arcs
     194      /// in a wraparound fashion and in the following minor iterations
     195      /// the best eligible arc is selected from this list.
     196      CANDIDATE_LIST,
     197
     198      /// The Altering Candidate List pivot rule.
     199      /// It is a modified version of the Candidate List method.
     200      /// It keeps only the several best eligible arcs from the former
     201      /// candidate list and extends this list in every iteration.
     202      ALTERING_LIST
     203    };
     204   
    186205  private:
    187206
    … …  
    221240    Node _psource, _ptarget;
    222241    Flow _pstflow;
    223     ProblemType _ptype;
     242    SupplyType _stype;
     243   
     244    Flow _sum_supply;
    224245
    225246    // Result maps
    … …  
    259280    int stem, par_stem, new_stem;
    260281    Flow delta;
     282
     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;
    261291
    262292  private:
    … …  
    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
    … …  
    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.
    … …  
    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);
    … …  
    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.
    … …  
    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);
    … …  
    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    ///
    … …  
    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);
    … …  
    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;
    … …  
    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.)
     815    ///
     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.
    822819    ///
    823820    /// \param s The source node.
    … …  
    837834    }
    838835   
    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
     836    /// \brief Set the type of the supply constraints.
     837    ///
     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    }
    … …  
    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
    … …  
    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
    … …  
    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    ///
    … …  
    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    ///
    … …  
    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
    … …  
    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;
    … …  
    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
    … …  
    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)
    … …  
    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.
    … …  
    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;
    … …  
    11131117          _node_id[n] = i;
    11141118          _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          _sum_supply += _supply[i];
     1120        }
     1121        valid_supply = (_stype == GEQ && _sum_supply <= 0) ||
     1122                       (_stype == LEQ && _sum_supply >= 0);
    11191123      } else {
    11201124        int i = 0;
    … …  
    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];
    1144         }
    1145         art_cost = (art_cost + 1) * _node_num;
    1146       }
    1147 
    1148       // Run Circulation to check if a feasible solution exists
    1149       typedef ConstMap<Arc, Flow> ConstArcMap;
    1150       ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
    1151       FlowNodeMap *csup = NULL;
    1152       bool local_csup = false;
    1153       if (_psupply) {
    1154         csup = _psupply;
    1155       } else {
    1156         csup = new FlowNodeMap(_graph, 0);
    1157         (*csup)[_psource] =  _pstflow;
    1158         (*csup)[_ptarget] = -_pstflow;
    1159         local_csup = true;
    1160       }
    1161       bool circ_result = false;
    1162       if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
    1163         // GEQ problem type
    1164         if (_plower) {
    1165           if (_pupper) {
    1166             Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
    1167               circ(_graph, *_plower, *_pupper, *csup);
    1168             circ_result = circ.run();
    1169           } else {
    1170             Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
    1171               circ(_graph, *_plower, inf_arc_map, *csup);
    1172             circ_result = circ.run();
    1173           }
    1174         } else {
    1175           if (_pupper) {
    1176             Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
    1177               circ(_graph, zero_arc_map, *_pupper, *csup);
    1178             circ_result = circ.run();
    1179           } else {
    1180             Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
    1181               circ(_graph, zero_arc_map, inf_arc_map, *csup);
    1182             circ_result = circ.run();
    1183           }
    1184         }
    1185       } else {
    1186         // LEQ problem type
    1187         typedef ReverseDigraph<const GR> RevGraph;
    1188         typedef NegMap<FlowNodeMap> NegNodeMap;
    1189         RevGraph rgraph(_graph);
    1190         NegNodeMap neg_csup(*csup);
    1191         if (_plower) {
    1192           if (_pupper) {
    1193             Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
    1194               circ(rgraph, *_plower, *_pupper, neg_csup);
    1195             circ_result = circ.run();
    1196           } else {
    1197             Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
    1198               circ(rgraph, *_plower, inf_arc_map, neg_csup);
    1199             circ_result = circ.run();
    1200           }
    1201         } else {
    1202           if (_pupper) {
    1203             Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
    1204               circ(rgraph, zero_arc_map, *_pupper, neg_csup);
    1205             circ_result = circ.run();
    1206           } else {
    1207             Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
    1208               circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
    1209             circ_result = circ.run();
    1210           }
    1211         }
    1212       }
    1213       if (local_csup) delete csup;
    1214       if (!circ_result) return false;
     1141          if (_cost[i] > ART_COST) ART_COST = _cost[i];
     1142        }
     1143        ART_COST = (ART_COST + 1) * _node_num;
     1144      }
    12151145
    12161146      // Set data for the artificial root node
    … …  
    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
    … …  
    12611191        } else {
    12621192          for (int i = 0; i != _arc_num; ++i)
    1263             _cap[i] = inf_cap;
     1193            _cap[i] = INF;
    12641194        }
    12651195        if (_pcost) {
    … …  
    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;
    … …  
    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      }
    … …  
    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;
    … …  
    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;
    … …  
    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) {
    … …  
    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
    … …  
    15521493        findJoinNode();
    15531494        bool change = findLeavingArc();
     1495        if (delta >= INF) return UNBOUNDED;
    15541496        changeFlow(change);
    15551497        if (change) {
    15561498          updateTreeStructure();
    15571499          updatePotential();
     1500        }
     1501      }
     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;
    15581517        }
    15591518      }
    … …  
    15751534      }
    15761535
    1577       return true;
     1536      return OPTIMAL;
    15781537    }
    15791538
Note: See TracChangeset for help on using the changeset viewer.