lemon-1.2: comparison lemon/network

equal deleted inserted replaced

-:88ad3fc9dbe1
+:d793a53e7068
 #include <limits>
 #include <algorithm>
 #include <lemon/core.h>
 #include <lemon/math.h>
-#include <lemon/maps.h>
-#include <lemon/circulation.h>
-#include <lemon/adaptors.h>
 namespace lemon {
 /// \addtogroup min_cost_flow
 /// @{
 /// simplex method directly for the minimum cost flow problem.
 /// It is one of the most efficient solution methods.
 ///
 /// In general this class is the fastest implementation available
 /// in LEMON for the minimum cost flow problem.
-/// Moreover it supports both direction of the supply/demand inequality
+/// Moreover it supports both directions of the supply/demand inequality
-/// constraints. For more information see \ref ProblemType.
+/// constraints. For more information see \ref SupplyType.
+///
+/// Most of the parameters of the problem (except for the digraph)
+/// can be given using separate functions, and the algorithm can be
+/// executed using the \ref run() function. If some parameters are not
+/// specified, then default values will be used.
 ///
 /// \tparam GR The digraph type the algorithm runs on.
 /// \tparam F The value type used for flow amounts, capacity bounds
 /// and supply values in the algorithm. By default it is \c int.
 /// \tparam C The value type used for costs and potentials in the
 typedef typename GR::template NodeMap<Cost> PotentialMap;
 #endif
 public:
-/// \brief Enum type for selecting the pivot rule.
+/// \brief Problem type constants for the \c run() function.
 ///
-/// Enum type for selecting the pivot rule for the \ref run()
+/// Enum type containing the problem type constants that can be
-/// function.
+/// returned by the \ref run() function of the algorithm.
-///
+enum ProblemType {
-/// \ref NetworkSimplex provides five different pivot rule
+/// The problem has no feasible solution (flow).
-/// implementations that significantly affect the running time
+INFEASIBLE,
-/// of the algorithm.
+/// The problem has optimal solution (i.e. it is feasible and
-/// By default \ref BLOCK_SEARCH "Block Search" is used, which
+/// bounded), and the algorithm has found optimal flow and node
-/// proved to be the most efficient and the most robust on various
+/// potentials (primal and dual solutions).
-/// test inputs according to our benchmark tests.
+OPTIMAL,
-/// However another pivot rule can be selected using the \ref run()
+/// The objective function of the problem is unbounded, i.e.
-/// function with the proper parameter.
+/// there is a directed cycle having negative total cost and
-enum PivotRule {
+/// infinite upper bound.
+UNBOUNDED
-/// The First Eligible pivot rule.
-/// The next eligible arc is selected in a wraparound fashion
-/// in every iteration.
-FIRST_ELIGIBLE,
-/// The Best Eligible pivot rule.
-/// The best eligible arc is selected in every iteration.
-BEST_ELIGIBLE,
-/// The Block Search pivot rule.
-/// A specified number of arcs are examined in every iteration
-/// in a wraparound fashion and the best eligible arc is selected
-/// from this block.
-BLOCK_SEARCH,
-/// The Candidate List pivot rule.
-/// In a major iteration a candidate list is built from eligible arcs
-/// in a wraparound fashion and in the following minor iterations
-/// the best eligible arc is selected from this list.
-CANDIDATE_LIST,
-/// The Altering Candidate List pivot rule.
-/// It is a modified version of the Candidate List method.
-/// It keeps only the several best eligible arcs from the former
-/// candidate list and extends this list in every iteration.
-ALTERING_LIST
 };
-/// \brief Enum type for selecting the problem type.
+/// \brief Constants for selecting the type of the supply constraints.
 ///
-/// Enum type for selecting the problem type, i.e. the direction of
+/// Enum type containing constants for selecting the supply type,
-/// the inequalities in the supply/demand constraints of the
+/// i.e. the direction of the inequalities in the supply/demand
-/// \ref min_cost_flow "minimum cost flow problem".
+/// constraints of the \ref min_cost_flow "minimum cost flow problem".
 ///
-/// The default problem type is \c GEQ, since this form is supported
+/// The default supply type is \c GEQ, since this form is supported
 /// by other minimum cost flow algorithms and the \ref Circulation
-/// algorithm as well.
+/// algorithm, as well.
-/// The \c LEQ problem type can be selected using the \ref problemType()
+/// The \c LEQ problem type can be selected using the \ref supplyType()
 /// function.
 ///
-/// Note that the equality form is a special case of both problem type.
+/// Note that the equality form is a special case of both supply types.
-enum ProblemType {
+enum SupplyType {
-/// This option means that there are "<em>greater or equal</em>"
+/// This option means that there are <em>"greater or equal"</em>
-/// constraints in the defintion, i.e. the exact formulation of the
+/// supply/demand constraints in the definition, i.e. the exact
-/// problem is the following.
+/// formulation of the problem is the following.
 /**
 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
 \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
 sup(u) \quad \forall u\in V \f]
 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
 /// satisfied.
 GEQ,
 /// It is just an alias for the \c GEQ option.
 CARRY_SUPPLIES = GEQ,
-/// This option means that there are "<em>less or equal</em>"
+/// This option means that there are <em>"less or equal"</em>
-/// constraints in the defintion, i.e. the exact formulation of the
+/// supply/demand constraints in the definition, i.e. the exact
-/// problem is the following.
+/// formulation of the problem is the following.
 /**
 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
 \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
 sup(u) \quad \forall u\in V \f]
 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
 /// nodes.
 LEQ,
 /// It is just an alias for the \c LEQ option.
 SATISFY_DEMANDS = LEQ
 };
+/// \brief Constants for selecting the pivot rule.
+///
+/// Enum type containing constants for selecting the pivot rule for
+/// the \ref run() function.
+///
+/// \ref NetworkSimplex provides five different pivot rule
+/// implementations that significantly affect the running time
+/// of the algorithm.
+/// By default \ref BLOCK_SEARCH "Block Search" is used, which
+/// proved to be the most efficient and the most robust on various
+/// test inputs according to our benchmark tests.
+/// However another pivot rule can be selected using the \ref run()
+/// function with the proper parameter.
+enum PivotRule {
+/// The First Eligible pivot rule.
+/// The next eligible arc is selected in a wraparound fashion
+/// in every iteration.
+FIRST_ELIGIBLE,
+/// The Best Eligible pivot rule.
+/// The best eligible arc is selected in every iteration.
+BEST_ELIGIBLE,
+/// The Block Search pivot rule.
+/// A specified number of arcs are examined in every iteration
+/// in a wraparound fashion and the best eligible arc is selected
+/// from this block.
+BLOCK_SEARCH,
+/// The Candidate List pivot rule.
+/// In a major iteration a candidate list is built from eligible arcs
+/// in a wraparound fashion and in the following minor iterations
+/// the best eligible arc is selected from this list.
+CANDIDATE_LIST,
+/// The Altering Candidate List pivot rule.
+/// It is a modified version of the Candidate List method.
+/// It keeps only the several best eligible arcs from the former
+/// candidate list and extends this list in every iteration.
+ALTERING_LIST
+};
 private:
 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 typedef typename GR::template ArcMap<Flow> FlowArcMap;
 CostArcMap *_pcost;
 FlowNodeMap *_psupply;
 bool _pstsup;
 Node _psource, _ptarget;
 Flow _pstflow;
-ProblemType _ptype;
+SupplyType _stype;
+Flow _sum_supply;
 // Result maps
 FlowMap *_flow_map;
 PotentialMap *_potential_map;
 bool _local_flow;
 // Temporary data used in the current pivot iteration
 int in_arc, join, u_in, v_in, u_out, v_out;
 int first, second, right, last;
 int stem, par_stem, new_stem;
 Flow delta;
+public:
+/// \brief Constant for infinite upper bounds (capacities).
+///
+/// Constant for infinite upper bounds (capacities).
+/// It is \c std::numeric_limits<Flow>::infinity() if available,
+/// \c std::numeric_limits<Flow>::max() otherwise.
+const Flow INF;
 private:
 // Implementation of the First Eligible pivot rule
 class FirstEligiblePivotRule
 ///
 /// \param graph The digraph the algorithm runs on.
 NetworkSimplex(const GR& graph) :
 _graph(graph),
 _plower(NULL), _pupper(NULL), _pcost(NULL),
-_psupply(NULL), _pstsup(false), _ptype(GEQ),
+_psupply(NULL), _pstsup(false), _stype(GEQ),
 _flow_map(NULL), _potential_map(NULL),
 _local_flow(false), _local_potential(false),
-_node_id(graph)
+_node_id(graph),
+INF(std::numeric_limits<Flow>::has_infinity ?
+std::numeric_limits<Flow>::infinity() :
+std::numeric_limits<Flow>::max())
 {
-LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
+// Check the value types
-std::numeric_limits<Flow>::is_signed,
+LEMON_ASSERT(std::numeric_limits<Flow>::is_signed,
-"The flow type of NetworkSimplex must be signed integer");
+"The flow type of NetworkSimplex must be signed");
-LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
+LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
-std::numeric_limits<Cost>::is_signed,
+"The cost type of NetworkSimplex must be signed");
-"The cost type of NetworkSimplex must be signed integer");
 }
 /// Destructor.
 ~NetworkSimplex() {
 if (_local_flow) delete _flow_map;
 /// @{
 /// \brief Set the lower bounds on the arcs.
 ///
 /// This function sets the lower bounds on the arcs.
-/// If neither this function nor \ref boundMaps() is used before
+/// If it is not used before calling \ref run(), the lower bounds
-/// calling \ref run(), the lower bounds will be set to zero
+/// will be set to zero on all arcs.
-/// on all arcs.
 ///
 /// \param map An arc map storing the lower bounds.
 /// Its \c Value type must be convertible to the \c Flow type
 /// of the algorithm.
 ///
 /// \return <tt>(*this)</tt>
-template <typename LOWER>
+template <typename LowerMap>
-NetworkSimplex& lowerMap(const LOWER& map) {
+NetworkSimplex& lowerMap(const LowerMap& map) {
 delete _plower;
 _plower = new FlowArcMap(_graph);
 for (ArcIt a(_graph); a != INVALID; ++a) {
 (*_plower)[a] = map[a];
 }
 }
 /// \brief Set the upper bounds (capacities) on the arcs.
 ///
 /// This function sets the upper bounds (capacities) on the arcs.
-/// If none of the functions \ref upperMap(), \ref capacityMap()
+/// If it is not used before calling \ref run(), the upper bounds
-/// and \ref boundMaps() is used before calling \ref run(),
+/// will be set to \ref INF on all arcs (i.e. the flow value will be
-/// the upper bounds (capacities) will be set to
+/// unbounded from above on each arc).
-/// \c std::numeric_limits<Flow>::max() on all arcs.
 ///
 /// \param map An arc map storing the upper bounds.
 /// Its \c Value type must be convertible to the \c Flow type
 /// of the algorithm.
 ///
 /// \return <tt>(*this)</tt>
-template<typename UPPER>
+template<typename UpperMap>
-NetworkSimplex& upperMap(const UPPER& map) {
+NetworkSimplex& upperMap(const UpperMap& map) {
 delete _pupper;
 _pupper = new FlowArcMap(_graph);
 for (ArcIt a(_graph); a != INVALID; ++a) {
 (*_pupper)[a] = map[a];
 }
 return *this;
 }
-/// \brief Set the upper bounds (capacities) on the arcs.
-///
-/// This function sets the upper bounds (capacities) on the arcs.
-/// It is just an alias for \ref upperMap().
-///
-/// \return <tt>(*this)</tt>
-template<typename CAP>
-NetworkSimplex& capacityMap(const CAP& map) {
-return upperMap(map);
-}
-/// \brief Set the lower and upper bounds on the arcs.
-///
-/// This function sets the lower and upper bounds on the arcs.
-/// If neither this function nor \ref lowerMap() is used before
-/// calling \ref run(), the lower bounds will be set to zero
-/// on all arcs.
-/// If none of the functions \ref upperMap(), \ref capacityMap()
-/// and \ref boundMaps() is used before calling \ref run(),
-/// the upper bounds (capacities) will be set to
-/// \c std::numeric_limits<Flow>::max() on all arcs.
-///
-/// \param lower An arc map storing the lower bounds.
-/// \param upper An arc map storing the upper bounds.
-///
-/// The \c Value type of the maps must be convertible to the
-/// \c Flow type of the algorithm.
-///
-/// \note This function is just a shortcut of calling \ref lowerMap()
-/// and \ref upperMap() separately.
-///
-/// \return <tt>(*this)</tt>
-template <typename LOWER, typename UPPER>
-NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
-return lowerMap(lower).upperMap(upper);
-}
 /// \brief Set the costs of the arcs.
 ///
 /// This function sets the costs of the arcs.
 /// If it is not used before calling \ref run(), the costs
 /// will be set to \c 1 on all arcs.
 /// \param map An arc map storing the costs.
 /// Its \c Value type must be convertible to the \c Cost type
 /// of the algorithm.
 ///
 /// \return <tt>(*this)</tt>
-template<typename COST>
+template<typename CostMap>
-NetworkSimplex& costMap(const COST& map) {
+NetworkSimplex& costMap(const CostMap& map) {
 delete _pcost;
 _pcost = new CostArcMap(_graph);
 for (ArcIt a(_graph); a != INVALID; ++a) {
 (*_pcost)[a] = map[a];
 }
 /// \param map A node map storing the supply values.
 /// Its \c Value type must be convertible to the \c Flow type
 /// of the algorithm.
 ///
 /// \return <tt>(*this)</tt>
-template<typename SUP>
+template<typename SupplyMap>
-NetworkSimplex& supplyMap(const SUP& map) {
+NetworkSimplex& supplyMap(const SupplyMap& map) {
 delete _psupply;
 _pstsup = false;
 _psupply = new FlowNodeMap(_graph);
 for (NodeIt n(_graph); n != INVALID; ++n) {
 (*_psupply)[n] = map[n];
 /// This function sets a single source node and a single target node
 /// and the required flow value.
 /// If neither this function nor \ref supplyMap() is used before
 /// calling \ref run(), the supply of each node will be set to zero.
 /// (It makes sense only if non-zero lower bounds are given.)
+///
+/// Using this function has the same effect as using \ref supplyMap()
+/// with such a map in which \c k is assigned to \c s, \c -k is
+/// assigned to \c t and all other nodes have zero supply value.
 ///
 /// \param s The source node.
 /// \param t The target node.
 /// \param k The required amount of flow from node \c s to node \c t
 /// (i.e. the supply of \c s and the demand of \c t).
 _ptarget = t;
 _pstflow = k;
 return *this;
 }
-/// \brief Set the problem type.
+/// \brief Set the type of the supply constraints.
 ///
-/// This function sets the problem type for the algorithm.
+/// This function sets the type of the supply/demand constraints.
-/// If it is not used before calling \ref run(), the \ref GEQ problem
+/// If it is not used before calling \ref run(), the \ref GEQ supply
 /// type will be used.
 ///
-/// For more information see \ref ProblemType.
+/// For more information see \ref SupplyType.
 ///
 /// \return <tt>(*this)</tt>
-NetworkSimplex& problemType(ProblemType problem_type) {
+NetworkSimplex& supplyType(SupplyType supply_type) {
-_ptype = problem_type;
+_stype = supply_type;
 return *this;
 }
 /// \brief Set the flow map.
 ///
 /// \brief Run the algorithm.
 ///
 /// This function runs the algorithm.
 /// The paramters can be specified using functions \ref lowerMap(),
-/// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
+/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
-/// \ref costMap(), \ref supplyMap(), \ref stSupply(),
+/// \ref supplyType(), \ref flowMap() and \ref potentialMap().
-/// \ref problemType(), \ref flowMap() and \ref potentialMap().
 /// For example,
 /// \code
 ///   NetworkSimplex<ListDigraph> ns(graph);
-///   ns.boundMaps(lower, upper).costMap(cost)
+///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
 ///     .supplyMap(sup).run();
 /// \endcode
 ///
 /// This function can be called more than once. All the parameters
 /// that have been given are kept for the next call, unless
 /// have to be set again. See \ref reset() for examples.
 ///
 /// \param pivot_rule The pivot rule that will be used during the
 /// algorithm. For more information see \ref PivotRule.
 ///
-/// \return \c true if a feasible flow can be found.
+/// \return \c INFEASIBLE if no feasible flow exists,
-bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
+/// \n \c OPTIMAL if the problem has optimal solution
-return init() && start(pivot_rule);
+/// (i.e. it is feasible and bounded), and the algorithm has found
+/// optimal flow and node potentials (primal and dual solutions),
+/// \n \c UNBOUNDED if the objective function of the problem is
+/// unbounded, i.e. there is a directed cycle having negative total
+/// cost and infinite upper bound.
+///
+/// \see ProblemType, PivotRule
+ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
+if (!init()) return INFEASIBLE;
+return start(pivot_rule);
 }
 /// \brief Reset all the parameters that have been given before.
 ///
 /// This function resets all the paramaters that have been given
 /// before using functions \ref lowerMap(), \ref upperMap(),
-/// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
+/// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(),
-/// \ref supplyMap(), \ref stSupply(), \ref problemType(),
 /// \ref flowMap() and \ref potentialMap().
 ///
 /// It is useful for multiple run() calls. If this function is not
 /// used, all the parameters given before are kept for the next
 /// \ref run() call.
 /// For example,
 /// \code
 ///   NetworkSimplex<ListDigraph> ns(graph);
 ///
 ///   // First run
-///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
+///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
 ///     .supplyMap(sup).run();
 ///
 ///   // Run again with modified cost map (reset() is not called,
 ///   // so only the cost map have to be set again)
 ///   cost[e] += 100;
 ///   ns.costMap(cost).run();
 ///
 ///   // Run again from scratch using reset()
 ///   // (the lower bounds will be set to zero on all arcs)
 ///   ns.reset();
-///   ns.capacityMap(cap).costMap(cost)
+///   ns.upperMap(capacity).costMap(cost)
 ///     .supplyMap(sup).run();
 /// \endcode
 ///
 /// \return <tt>(*this)</tt>
 NetworkSimplex& reset() {
 _plower = NULL;
 _pupper = NULL;
 _pcost = NULL;
 _psupply = NULL;
 _pstsup = false;
-_ptype = GEQ;
+_stype = GEQ;
 if (_local_flow) delete _flow_map;
 if (_local_potential) delete _potential_map;
 _flow_map = NULL;
 _potential_map = NULL;
 _local_flow = false;
 /// @{
 /// \brief Return the total cost of the found flow.
 ///
 /// This function returns the total cost of the found flow.
-/// The complexity of the function is O(e).
+/// Its complexity is O(e).
 ///
 /// \note The return type of the function can be specified as a
 /// template parameter. For example,
 /// \code
 ///   ns.totalCost<double>();
 /// It is useful if the total cost cannot be stored in the \c Cost
 /// type of the algorithm, which is the default return type of the
 /// function.
 ///
 /// \pre \ref run() must be called before using this function.
-template <typename Num>
+template <typename Value>
-Num totalCost() const {
+Value totalCost() const {
-Num c = 0;
+Value c = 0;
 if (_pcost) {
 for (ArcIt e(_graph); e != INVALID; ++e)
 c += (*_flow_map)[e] * (*_pcost)[e];
 } else {
 for (ArcIt e(_graph); e != INVALID; ++e)
 /// \brief Return a const reference to the potential map
 /// (the dual solution).
 ///
 /// This function returns a const reference to a node map storing
 /// the found potentials, which form the dual solution of the
-/// \ref min_cost_flow "minimum cost flow" problem.
+/// \ref min_cost_flow "minimum cost flow problem".
 ///
 /// \pre \ref run() must be called before using this function.
 const PotentialMap& potentialMap() const {
 return *_potential_map;
 }
 _last_succ.resize(all_node_num);
 _state.resize(all_arc_num);
 // Initialize node related data
 bool valid_supply = true;
-Flow sum_supply = 0;
+_sum_supply = 0;
 if (!_pstsup && !_psupply) {
 _pstsup = true;
 _psource = _ptarget = NodeIt(_graph);
 _pstflow = 0;
 }
 if (_psupply) {
 int i = 0;
 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
 _node_id[n] = i;
 _supply[i] = (*_psupply)[n];
-sum_supply += _supply[i];
+_sum_supply += _supply[i];
 }
-valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
+valid_supply = (_stype == GEQ && _sum_supply <= 0) ||
-(_ptype == LEQ && sum_supply >= 0);
+(_stype == LEQ && _sum_supply >= 0);
 } else {
 int i = 0;
 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
 _node_id[n] = i;
 _supply[i] = 0;
 _supply[_node_id[_psource]] =  _pstflow;
 _supply[_node_id[_ptarget]] = -_pstflow;
 }
 if (!valid_supply) return false;
-// Infinite capacity value
-Flow inf_cap =
-std::numeric_limits<Flow>::has_infinity ?
-std::numeric_limits<Flow>::infinity() :
-std::numeric_limits<Flow>::max();
 // Initialize artifical cost
-Cost art_cost;
+Cost ART_COST;
 if (std::numeric_limits<Cost>::is_exact) {
-art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
+ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
 } else {
-art_cost = std::numeric_limits<Cost>::min();
+ART_COST = std::numeric_limits<Cost>::min();
 for (int i = 0; i != _arc_num; ++i) {
-if (_cost[i] > art_cost) art_cost = _cost[i];
+if (_cost[i] > ART_COST) ART_COST = _cost[i];
 }
-art_cost = (art_cost + 1) * _node_num;
+ART_COST = (ART_COST + 1) * _node_num;
 }
-// Run Circulation to check if a feasible solution exists
-typedef ConstMap<Arc, Flow> ConstArcMap;
-ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
-FlowNodeMap *csup = NULL;
-bool local_csup = false;
-if (_psupply) {
-csup = _psupply;
-} else {
-csup = new FlowNodeMap(_graph, 0);
-(*csup)[_psource] =  _pstflow;
-(*csup)[_ptarget] = -_pstflow;
-local_csup = true;
-}
-bool circ_result = false;
-if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
-// GEQ problem type
-if (_plower) {
-if (_pupper) {
-Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
-circ(_graph, *_plower, *_pupper, *csup);
-circ_result = circ.run();
-} else {
-Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
-circ(_graph, *_plower, inf_arc_map, *csup);
-circ_result = circ.run();
-}
-} else {
-if (_pupper) {
-Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
-circ(_graph, zero_arc_map, *_pupper, *csup);
-circ_result = circ.run();
-} else {
-Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
-circ(_graph, zero_arc_map, inf_arc_map, *csup);
-circ_result = circ.run();
-}
-}
-} else {
-// LEQ problem type
-typedef ReverseDigraph<const GR> RevGraph;
-typedef NegMap<FlowNodeMap> NegNodeMap;
-RevGraph rgraph(_graph);
-NegNodeMap neg_csup(*csup);
-if (_plower) {
-if (_pupper) {
-Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
-circ(rgraph, *_plower, *_pupper, neg_csup);
-circ_result = circ.run();
-} else {
-Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
-circ(rgraph, *_plower, inf_arc_map, neg_csup);
-circ_result = circ.run();
-}
-} else {
-if (_pupper) {
-Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
-circ(rgraph, zero_arc_map, *_pupper, neg_csup);
-circ_result = circ.run();
-} else {
-Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
-circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
-circ_result = circ.run();
-}
-}
-}
-if (local_csup) delete csup;
-if (!circ_result) return false;
 // Set data for the artificial root node
 _root = _node_num;
 _parent[_root] = -1;
 _pred[_root] = -1;
 _thread[_root] = 0;
 _rev_thread[0] = _root;
 _succ_num[_root] = all_node_num;
 _last_succ[_root] = _root - 1;
-_supply[_root] = -sum_supply;
+_supply[_root] = -_sum_supply;
-if (sum_supply < 0) {
+if (_sum_supply < 0) {
-_pi[_root] = -art_cost;
+_pi[_root] = -ART_COST;
 } else {
-_pi[_root] = art_cost;
+_pi[_root] = ART_COST;
 }
 // Store the arcs in a mixed order
 int k = std::max(int(std::sqrt(double(_arc_num))), 10);
 int i = 0;
 if (_pupper) {
 for (int i = 0; i != _arc_num; ++i)
 _cap[i] = (*_pupper)[_arc_ref[i]];
 } else {
 for (int i = 0; i != _arc_num; ++i)
-_cap[i] = inf_cap;
+_cap[i] = INF;
 }
 if (_pcost) {
 for (int i = 0; i != _arc_num; ++i)
 _cost[i] = (*_pcost)[_arc_ref[i]];
 } else {
 // Remove non-zero lower bounds
 if (_plower) {
 for (int i = 0; i != _arc_num; ++i) {
 Flow c = (*_plower)[_arc_ref[i]];
-if (c != 0) {
+if (c > 0) {
-_cap[i] -= c;
+if (_cap[i] < INF) _cap[i] -= c;
+_supply[_source[i]] -= c;
+_supply[_target[i]] += c;
+}
+else if (c < 0) {
+if (_cap[i] < INF + c) {
+_cap[i] -= c;
+} else {
+_cap[i] = INF;
+}
 _supply[_source[i]] -= c;
 _supply[_target[i]] += c;
 }
 }
 }
 _rev_thread[u + 1] = u;
 _succ_num[u] = 1;
 _last_succ[u] = u;
 _parent[u] = _root;
 _pred[u] = e;
-_cost[e] = art_cost;
+_cost[e] = ART_COST;
-_cap[e] = inf_cap;
+_cap[e] = INF;
 _state[e] = STATE_TREE;
-if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
+if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
 _flow[e] = _supply[u];
 _forward[u] = true;
-_pi[u] = -art_cost + _pi[_root];
+_pi[u] = -ART_COST + _pi[_root];
 } else {
 _flow[e] = -_supply[u];
 _forward[u] = false;
-_pi[u] = art_cost + _pi[_root];
+_pi[u] = ART_COST + _pi[_root];
 }
 }
 return true;
 }
 int e;
 // Search the cycle along the path form the first node to the root
 for (int u = first; u != join; u = _parent[u]) {
 e = _pred[u];
-d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
+d = _forward[u] ?
+_flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
 if (d < delta) {
 delta = d;
 u_out = u;
 result = 1;
 }
 }
 // Search the cycle along the path form the second node to the root
 for (int u = second; u != join; u = _parent[u]) {
 e = _pred[u];
-d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
+d = _forward[u] ?
+(_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
 if (d <= delta) {
 delta = d;
 u_out = u;
 result = 2;
 }
 _pi[u] += sigma;
 }
 }
 // Execute the algorithm
-bool start(PivotRule pivot_rule) {
+ProblemType start(PivotRule pivot_rule) {
 // Select the pivot rule implementation
 switch (pivot_rule) {
 case FIRST_ELIGIBLE:
 return start<FirstEligiblePivotRule>();
 case BEST_ELIGIBLE:
 case CANDIDATE_LIST:
 return start<CandidateListPivotRule>();
 case ALTERING_LIST:
 return start<AlteringListPivotRule>();
 }
-return false;
+return INFEASIBLE; // avoid warning
 }
 template <typename PivotRuleImpl>
-bool start() {
+ProblemType start() {
 PivotRuleImpl pivot(*this);
 // Execute the Network Simplex algorithm
 while (pivot.findEnteringArc()) {
 findJoinNode();
 bool change = findLeavingArc();
+if (delta >= INF) return UNBOUNDED;
 changeFlow(change);
 if (change) {
 updateTreeStructure();
 updatePotential();
+}
+}
+// Check feasibility
+if (_sum_supply < 0) {
+for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
+}
+}
+else if (_sum_supply > 0) {
+for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
+}
+}
+else {
+for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+if (_flow[e] != 0) return INFEASIBLE;
 }
 }
 // Copy flow values to _flow_map
 if (_plower) {
 // Copy potential values to _potential_map
 for (NodeIt n(_graph); n != INVALID; ++n) {
 _potential_map->set(n, _pi[_node_id[n]]);
 }
-return true;
+return OPTIMAL;
 }
 }; //class NetworkSimplex
 ///@}

changeset 640	6c408d864fa1
parent 618	b95898314e09
child 641	756a5ec551c8