LEMON/LEMON-main Changeset - r640:6c408d864fa1

Repositories » LEMON » LEMON-main

Summary
Changelog
Files
Switch To
- loading...
Options
- Lightweight changelog
- Search
Pull Requests

Changeset - r640:6c408d864fa1

« 626:58357e

641:756a5e »

r640:6c408d864fa1

Wed, 29 Apr 2009 03:15:24

[Not Reviewed]

0 comments (0 inline)

kpeter (Peter Kovacs)
kpeter@inf.elte.hu

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.

0 4 0

default

4 files changed with 350 insertions and 332 deletions:

doc/groups.dox

lemon/network_simplex.h

203

244

test/min_cost_flow_test.cc

135

tools/dimacs-solver.cc

↑ Collapse diff ↑

doc/groups.dox

Show inline comments

@@ -347,27 +347,27 @@
 		This group contains the algorithms for finding minimum cost flows and
 		circulations.
 		The \e minimum \e cost \e flow \e problem is to find a feasible flow of
 		minimum total cost from a set of supply nodes to a set of demand nodes
 		in a network with capacity constraints (lower and upper bounds)
 		and arc costs.
 		Formally, let \f$G=(V,A)\f$ be a digraph,
 		\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
 		Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
 		\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
 		upper bounds for the flow values on the arcs, for which
 		\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
 		\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
-		on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
+		\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
 		\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
 		on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
 		signed supply values of the nodes.
 		If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
 		supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
 		\f$-sup(u)\f$ demand.
-		A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
 		A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
 		of the following optimization problem.
 		\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]
 		The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
@@ -399,34 +399,34 @@
 		    sup(u) \quad \forall u\in V. \f]
 		However if the sum of the supply values is zero, then these two problems
 		are equivalent. So if you need the equality form, you have to ensure this
 		additional contraint for the algorithms.
 		The dual solution of the minimum cost flow problem is represented by node
 		potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
-		An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
 		An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
 		is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
 		node potentials the following \e complementary \e slackness optimality
 		conditions hold.
 		 - For all \f$uv\in A\f$ arcs:
 		   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
 		   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
 		   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
 		 - For all \f$u\in V\f$:
+		 - For all \f$u\in V\f$ nodes:
 		   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
 		     then \f$\pi(u)=0\f$.
 		Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
-		\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
+		\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
 		\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
 		All algorithms provide dual solution (node potentials) as well
+		All algorithms provide dual solution (node potentials) as well,
 		if an optimal flow is found.
 		LEMON contains several algorithms for solving minimum cost flow problems.
 		 - \ref NetworkSimplex Primal Network Simplex algorithm with various
 		   pivot strategies.
 		 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
 		   cost scaling.
 		 - \ref CapacityScaling Successive Shortest %Path algorithm with optional

lemon/network_simplex.h

Show inline comments

@@ -25,19 +25,16 @@
 		/// \brief Network Simplex algorithm for finding a minimum cost flow.
 		#include <vector>
 		#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
 		  /// @{
 		  /// \brief Implementation of the primal Network Simplex algorithm
 		  /// for finding a \ref min_cost_flow "minimum cost flow".
@@ -45,18 +42,23 @@
 		  /// \ref NetworkSimplex implements the primal Network Simplex algorithm
 		  /// for finding a \ref min_cost_flow "minimum cost flow".
 		  /// This algorithm is a specialized version of the linear programming
 		  /// 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
 		  /// constraints. For more information see \ref ProblemType.
 		  /// Moreover it supports both directions of the supply/demand inequality
 		  /// 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
 		  /// algorithm. By default it is the same as \c F.
 		  ///
 		  /// \warning Both value types must be signed and all input data must
@@ -83,21 +85,90 @@
 		    /// The type of the flow map
 		    typedef typename GR::template ArcMap<Flow> FlowMap;
 		    /// The type of the potential map
 		    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
 		    /// returned by the \ref run() function of the algorithm.
 		    enum ProblemType {
 		      /// The problem has no feasible solution (flow).
 		      INFEASIBLE,
 		      /// The problem has optimal solution (i.e. it is feasible and
 		      /// bounded), and the algorithm has found optimal flow and node
 		      /// potentials (primal and dual solutions).
 		      OPTIMAL,
 		      /// The objective function of the problem is unbounded, i.e.
 		      /// there is a directed cycle having negative total cost and
 		      /// infinite upper bound.
 		      UNBOUNDED
 		    };
 		    /// \brief Constants for selecting the type of the supply constraints.
 		    ///
 		    /// Enum type containing constants for selecting the supply type,
 		    /// i.e. the direction of the inequalities in the supply/demand
 		    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
 		    ///
 		    /// 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.
 		    /// The \c LEQ problem type can be selected using the \ref supplyType()
 		    /// function.
 		    ///
 		    /// Note that the equality form is a special case of both supply types.
 		    enum SupplyType {
 		      /// This option means that there are <em>"greater or equal"</em>
 		      /// supply/demand constraints in the definition, i.e. the exact
 		      /// 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]
 		      */
 		      /// It means that the total demand must be greater or equal to the
 		      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
 		      /// negative) and all the supplies have to be carried out from
 		      /// the supply nodes, but there could be demands that are not
 		      /// 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>
 		      /// supply/demand constraints in the definition, i.e. the exact
 		      /// 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]
 		      */
 		      /// It means that the total demand must be less or equal to the
 		      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
 		      /// positive) and all the demands have to be satisfied, but there
 		      /// could be supplies that are not carried out from the supply
 		      /// 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.
@@ -126,68 +197,16 @@
 		      /// 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.
 		    ///
 		    /// Enum type for selecting the problem type, i.e. the direction of
 		    /// the inequalities in the supply/demand constraints of the
 		    /// \ref min_cost_flow "minimum cost flow problem".
 		    ///
 		    /// The default problem type is \c GEQ, since this form is supported
 		    /// by other minimum cost flow algorithms and the \ref Circulation
 		    /// algorithm as well.
 		    /// The \c LEQ problem type can be selected using the \ref problemType()
 		    /// function.
 		    ///
 		    /// Note that the equality form is a special case of both problem type.
 		    enum ProblemType {
 		      /// This option means that there are "<em>greater or equal</em>"
 		      /// constraints in the defintion, i.e. the exact 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]
 		      */
 		      /// It means that the total demand must be greater or equal to the
 		      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
 		      /// negative) and all the supplies have to be carried out from
 		      /// the supply nodes, but there could be demands that are not
 		      /// 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>"
 		      /// constraints in the defintion, i.e. the exact 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]
 		      */
 		      /// It means that the total demand must be less or equal to the
 		      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
 		      /// positive) and all the demands have to be satisfied, but there
 		      /// could be supplies that are not carried out from the supply
 		      /// nodes.
 		      LEQ,
 		      /// It is just an alias for the \c LEQ option.
 		      SATISFY_DEMANDS = LEQ
 		    };
 		  private:
 		    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 		    typedef typename GR::template ArcMap<Flow> FlowArcMap;
 		    typedef typename GR::template ArcMap<Cost> CostArcMap;
 		    typedef typename GR::template NodeMap<Flow> FlowNodeMap;
@@ -215,17 +234,19 @@
 		    // Parameters of the problem
 		    FlowArcMap *_plower;
 		    FlowArcMap *_pupper;
 		    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;
 		    bool _local_potential;
 		    // Data structures for storing the digraph
@@ -254,16 +275,25 @@
 		    int _root;
 		    // 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
+		    {
 		    private:
 		      // References to the NetworkSimplex class
@@ -656,27 +686,29 @@
 		    /// \brief Constructor.
 		    ///
 		    /// The constructor of the class.
 		    ///
 		    /// \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 &&
 		                   std::numeric_limits<Flow>::is_signed,
 		        "The flow type of NetworkSimplex must be signed integer");
 		      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
 		                   std::numeric_limits<Cost>::is_signed,
 		        "The cost type of NetworkSimplex must be signed integer");
 		      // Check the value types
 		      LEMON_ASSERT(std::numeric_limits<Flow>::is_signed,
 		        "The flow type of NetworkSimplex must be signed");
 		      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
 		        "The cost type of NetworkSimplex must be signed");
+		    }
 		    /// Destructor.
 		    ~NetworkSimplex() {
 		      if (_local_flow) delete _flow_map;
 		      if (_local_potential) delete _potential_map;
+		    }
@@ -684,108 +716,69 @@
 		    /// The parameters of the algorithm can be specified using these
 		    /// functions.
 		    /// @{
 		    /// \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
 		    /// calling \ref run(), the lower bounds will be set to zero
-		    /// on all arcs.
+		    /// If it is not used before calling \ref run(), the lower bounds
 		    /// will be set to zero 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>
 		    NetworkSimplex& lowerMap(const LOWER& map) {
 		    template <typename LowerMap>
 		    NetworkSimplex& lowerMap(const LowerMap& map) {
 		      delete _plower;
 		      _plower = new FlowArcMap(_graph);
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        (*_plower)[a] = map[a];
+		      }
 		      return *this;
+		    }
 		    /// \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()
 		    /// 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.
 		    /// If it is not used before calling \ref run(), the upper bounds
 		    /// will be set to \ref INF on all arcs (i.e. the flow value will be
 		    /// unbounded from above on each arc).
 		    ///
 		    /// \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>
 		    NetworkSimplex& upperMap(const UPPER& map) {
 		    template<typename UpperMap>
 		    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>
 		    NetworkSimplex& costMap(const COST& map) {
 		    template<typename CostMap>
 		    NetworkSimplex& costMap(const CostMap& map) {
 		      delete _pcost;
 		      _pcost = new CostArcMap(_graph);
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        (*_pcost)[a] = map[a];
+		      }
 		      return *this;
+		    }
@@ -796,18 +789,18 @@
 		    /// calling \ref run(), the supply of each node will be set to zero.
 		    /// (It makes sense only if non-zero lower bounds are given.)
 		    ///
 		    /// \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>
 		    NetworkSimplex& supplyMap(const SUP& map) {
 		    template<typename SupplyMap>
 		    NetworkSimplex& supplyMap(const SupplyMap& map) {
 		      delete _psupply;
 		      _pstsup = false;
 		      _psupply = new FlowNodeMap(_graph);
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        (*_psupply)[n] = map[n];
+		      }
 		      return *this;
+		    }
@@ -815,43 +808,47 @@
 		    /// \brief Set single source and target nodes and a supply value.
 		    ///
 		    /// 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).
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
 		      delete _psupply;
 		      _psupply = NULL;
 		      _pstsup = true;
 		      _psource = s;
 		      _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.
 		    /// If it is not used before calling \ref run(), the \ref GEQ problem
 		    /// This function sets the type of the supply/demand constraints.
 		    /// 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) {
 		      _ptype = problem_type;
 		    NetworkSimplex& supplyType(SupplyType supply_type) {
 		      _stype = supply_type;
 		      return *this;
+		    }
 		    /// \brief Set the flow map.
 		    ///
 		    /// This function sets the flow map.
 		    /// If it is not used before calling \ref run(), an instance will
 		    /// be allocated automatically. The destructor deallocates this
@@ -891,83 +888,90 @@
 		    /// The algorithm can be executed using \ref run().
 		    /// @{
 		    /// \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 costMap(), \ref supplyMap(), \ref stSupply(),
-		    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
+		    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
 		    /// \ref supplyType(), \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
 		    /// \ref reset() is called, thus only the modified parameters
 		    /// 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.
 		    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
-		      return init() && start(pivot_rule);
+		    /// \return \c INFEASIBLE if no feasible flow exists,
 		    /// \n \c OPTIMAL if the problem has optimal solution
 		    /// (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 supplyMap(), \ref stSupply(), \ref problemType(),
 		    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(),
 		    /// \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() {
 		      delete _plower;
 		      delete _pupper;
 		      delete _pcost;
 		      delete _psupply;
 		      _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;
 		      _local_potential = false;
 		      return *this;
@@ -980,31 +984,31 @@
 		    /// functions.\n
 		    /// The \ref run() function must be called before using them.
 		    /// @{
 		    /// \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>();
 		    /// \endcode
 		    /// 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>
 		    Num totalCost() const {
-		      Num c = 0;
+		    template <typename Value>
 		    Value totalCost() const {
 		      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)
 		          c += (*_flow_map)[e];
+		      }
 		      return c;
@@ -1045,17 +1049,17 @@
 		      return (*_potential_map)[n];
+		    }
 		    /// \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;
+		    }
 		    /// @}
@@ -1096,141 +1100,67 @@
 		      _thread.resize(all_node_num);
 		      _rev_thread.resize(all_node_num);
 		      _succ_num.resize(all_node_num);
 		      _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) ||
 		                       (_ptype == LEQ && sum_supply >= 0);
 		        valid_supply = (_stype == GEQ && _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;
 		      if (sum_supply < 0) {
-		        _pi[_root] = -art_cost;
+		      _supply[_root] = -_sum_supply;
 		      if (_sum_supply < 0) {
 		        _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;
 		      for (ArcIt e(_graph); e != INVALID; ++e) {
 		        _arc_ref[i] = e;
 		        if ((i += k) >= _arc_num) i = (i % k) + 1;
@@ -1255,58 +1185,67 @@
 		          _flow[i] = 0;
 		          _state[i] = STATE_LOWER;
+		        }
 		        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 {
 		          for (int i = 0; i != _arc_num; ++i)
 		            _cost[i] = 1;
+		        }
+		      }
 		      // Remove non-zero lower bounds
 		      if (_plower) {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Flow c = (*_plower)[_arc_ref[i]];
 		          if (c != 0) {
 		            _cap[i] -= c;
 		          if (c > 0) {
 		            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;
+		          }
+		        }
+		      }
 		      // Add artificial arcs and initialize the spanning tree data structure
 		      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 		        _thread[u] = u + 1;
 		        _rev_thread[u + 1] = u;
 		        _succ_num[u] = 1;
 		        _last_succ[u] = u;
 		        _parent[u] = _root;
 		        _pred[u] = e;
 		        _cost[e] = art_cost;
 		        _cap[e] = inf_cap;
 		        _cost[e] = ART_COST;
 		        _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;
+		    }
 		    // Find the join node
 		    void findJoinNode() {
@@ -1337,27 +1276,29 @@
 		      delta = _cap[in_arc];
 		      int result = 0;
 		      Flow d;
 		      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;
+		        }
+		      }
 		      if (result == 1) {
@@ -1521,47 +1462,65 @@
 		      // Update potentials in the subtree, which has been moved
 		      int end = _thread[_last_succ[u_in]];
 		      for (int u = u_in; u != end; u = _thread[u]) {
 		        _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:
 		          return start<BestEligiblePivotRule>();
 		        case BLOCK_SEARCH:
 		          return start<BlockSearchPivotRule>();
 		        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) {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Arc e = _arc_ref[i];
 		          _flow_map->set(e, (*_plower)[e] + _flow[i]);
+		        }
 		      } else {
@@ -1569,17 +1528,17 @@
 		          _flow_map->set(_arc_ref[i], _flow[i]);
+		        }
+		      }
 		      // 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
 		  ///@}
 		} //namespace lemon

test/min_cost_flow_test.cc

Show inline comments

@@ -13,75 +13,76 @@
 		 * This software is provided "AS IS" with no warranty of any kind,
 		 * express or implied, and with no claim as to its suitability for any
 		 * purpose.
+		 *
 		 */
 		#include <iostream>
 		#include <fstream>
 		#include <limits>
 		#include <lemon/list_graph.h>
 		#include <lemon/lgf_reader.h>
 		#include <lemon/network_simplex.h>
 		#include <lemon/concepts/digraph.h>
 		#include <lemon/concept_check.h>
 		#include "test_tools.h"
 		using namespace lemon;
 		char test_lgf[] =
 		  "@nodes\n"
 		  "label  sup1 sup2 sup3 sup4 sup5\n"
 		  "    1    20   27    0   20   30\n"
 		  "    2    -4    0    0   -8   -3\n"
 		  "    3     0    0    0    0    0\n"
 		  "    4     0    0    0    0    0\n"
 		  "    5     9    0    0    6   11\n"
 		  "    6    -6    0    0   -5   -6\n"
 		  "    7     0    0    0    0    0\n"
 		  "    8     0    0    0    0    3\n"
 		  "    9     3    0    0    0    0\n"
 		  "   10    -2    0    0   -7   -2\n"
 		  "   11     0    0    0  -10    0\n"
 		  "   12   -20  -27    0  -30  -20\n"
 		  "\n"
 		  "label  sup1 sup2 sup3 sup4 sup5 sup6\n"
 		  "    1    20   27    0   30   20   30\n"
 		  "    2    -4    0    0    0   -8   -3\n"
 		  "    3     0    0    0    0    0    0\n"
 		  "    4     0    0    0    0    0    0\n"
 		  "    5     9    0    0    0    6   11\n"
 		  "    6    -6    0    0    0   -5   -6\n"
 		  "    7     0    0    0    0    0    0\n"
 		  "    8     0    0    0    0    0    3\n"
 		  "    9     3    0    0    0    0    0\n"
 		  "   10    -2    0    0    0   -7   -2\n"
 		  "   11     0    0    0    0  -10    0\n"
 		  "   12   -20  -27    0  -30  -30  -20\n"
 		  "\n"
 		  "@arcs\n"
 		  "       cost  cap low1 low2\n"
 		  " 1  2    70   11    0    8\n"
 		  " 1  3   150    3    0    1\n"
 		  " 1  4    80   15    0    2\n"
 		  " 2  8    80   12    0    0\n"
 		  " 3  5   140    5    0    3\n"
 		  " 4  6    60   10    0    1\n"
 		  " 4  7    80    2    0    0\n"
 		  " 4  8   110    3    0    0\n"
 		  " 5  7    60   14    0    0\n"
 		  " 5 11   120   12    0    0\n"
 		  " 6  3     0    3    0    0\n"
 		  " 6  9   140    4    0    0\n"
 		  " 6 10    90    8    0    0\n"
 		  " 7  1    30    5    0    0\n"
 		  " 8 12    60   16    0    4\n"
 		  " 9 12    50    6    0    0\n"
 		  "10 12    70   13    0    5\n"
 		  "10  2   100    7    0    0\n"
 		  "10  7    60   10    0    0\n"
 		  "11 10    20   14    0    6\n"
 		  "12 11    30   10    0    0\n"
 		  "       cost  cap low1 low2 low3\n"
 		  " 1  2    70   11    0    8    8\n"
 		  " 1  3   150    3    0    1    0\n"
 		  " 1  4    80   15    0    2    2\n"
 		  " 2  8    80   12    0    0    0\n"
 		  " 3  5   140    5    0    3    1\n"
 		  " 4  6    60   10    0    1    0\n"
 		  " 4  7    80    2    0    0    0\n"
 		  " 4  8   110    3    0    0    0\n"
 		  " 5  7    60   14    0    0    0\n"
 		  " 5 11   120   12    0    0    0\n"
 		  " 6  3     0    3    0    0    0\n"
 		  " 6  9   140    4    0    0    0\n"
 		  " 6 10    90    8    0    0    0\n"
 		  " 7  1    30    5    0    0   -5\n"
 		  " 8 12    60   16    0    4    3\n"
 		  " 9 12    50    6    0    0    0\n"
 		  "10 12    70   13    0    5    2\n"
 		  "10  2   100    7    0    0    0\n"
 		  "10  7    60   10    0    0   -3\n"
 		  "11 10    20   14    0    6  -20\n"
 		  "12 11    30   10    0    0  -10\n"
 		  "\n"
 		  "@attributes\n"
 		  "source 1\n"
 		  "target 12\n";
-		enum ProblemType {
+		enum SupplyType {
 		  EQ,
 		  GEQ,
 		  LEQ
 		};
 		// Check the interface of an MCF algorithm
 		template <typename GR, typename Flow, typename Cost>
 		class McfClassConcept
@@ -93,34 +94,34 @@
 		    void constraints() {
 		      checkConcept<concepts::Digraph, GR>();
 		      MCF mcf(g);
 		      b = mcf.reset()
 		             .lowerMap(lower)
 		             .upperMap(upper)
 		             .capacityMap(upper)
 		             .boundMaps(lower, upper)
 		             .costMap(cost)
 		             .supplyMap(sup)
 		             .stSupply(n, n, k)
 		             .flowMap(flow)
 		             .potentialMap(pot)
 		             .run();
 		      const MCF& const_mcf = mcf;
 		      const typename MCF::FlowMap &fm = const_mcf.flowMap();
 		      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
-		      v = const_mcf.totalCost();
+		      c = const_mcf.totalCost();
 		      double x = const_mcf.template totalCost<double>();
 		      v = const_mcf.flow(a);
-		      v = const_mcf.potential(n);
+		      c = const_mcf.potential(n);
 		      v = const_mcf.INF;
 		      ignore_unused_variable_warning(fm);
 		      ignore_unused_variable_warning(pm);
 		      ignore_unused_variable_warning(x);
+		    }
 		    typedef typename GR::Node Node;
 		    typedef typename GR::Arc Arc;
@@ -132,31 +133,32 @@
 		    const FAM &lower;
 		    const FAM &upper;
 		    const CAM &cost;
 		    const NM &sup;
 		    const Node &n;
 		    const Arc &a;
 		    const Flow &k;
 		    Flow v;
 		    Cost c;
 		    bool b;
 		    typename MCF::FlowMap &flow;
 		    typename MCF::PotentialMap &pot;
 		  };
 		};
 		// Check the feasibility of the given flow (primal soluiton)
 		template < typename GR, typename LM, typename UM,
 		           typename SM, typename FM >
 		bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
 		                const SM& supply, const FM& flow,
-		                ProblemType type = EQ )
+		                SupplyType type = EQ )
+		{
 		  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 		  for (ArcIt e(gr); e != INVALID; ++e) {
 		    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
+		  }
 		  for (NodeIt n(gr); n != INVALID; ++n) {
@@ -203,26 +205,27 @@
+		  }
 		  return opt;
+		}
 		// Run a minimum cost flow algorithm and check the results
 		template < typename MCF, typename GR,
 		           typename LM, typename UM,
 		           typename CM, typename SM >
 		void checkMcf( const MCF& mcf, bool mcf_result,
 		           typename CM, typename SM,
 		           typename PT >
 		void checkMcf( const MCF& mcf, PT mcf_result,
 		               const GR& gr, const LM& lower, const UM& upper,
 		               const CM& cost, const SM& supply,
-		               bool result, typename CM::Value total,
+		               PT result, bool optimal, typename CM::Value total,
 		               const std::string &test_id = "",
-		               ProblemType type = EQ )
+		               SupplyType type = EQ )
+		{
 		  check(mcf_result == result, "Wrong result " + test_id);
-		  if (result) {
+		  if (optimal) {
 		    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
 		          "The flow is not feasible " + test_id);
 		    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
 		    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
 		                         mcf.potentialMap()),
 		          "Wrong potentials " + test_id);
+		  }
+		}
@@ -239,96 +242,152 @@
+		  }
 		  // Run various MCF tests
 		  typedef ListDigraph Digraph;
 		  DIGRAPH_TYPEDEFS(ListDigraph);
 		  // Read the test digraph
 		  Digraph gr;
 		  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
 		  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
 		  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
 		  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
 		  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
 		  Node v, w;
 		  std::istringstream input(test_lgf);
 		  DigraphReader<Digraph>(gr, input)
 		    .arcMap("cost", c)
 		    .arcMap("cap", u)
 		    .arcMap("low1", l1)
 		    .arcMap("low2", l2)
 		    .arcMap("low3", l3)
 		    .nodeMap("sup1", s1)
 		    .nodeMap("sup2", s2)
 		    .nodeMap("sup3", s3)
 		    .nodeMap("sup4", s4)
 		    .nodeMap("sup5", s5)
 		    .nodeMap("sup6", s6)
 		    .node("source", v)
 		    .node("target", w)
 		    .run();
 		  // Build a test digraph for testing negative costs
 		  Digraph ngr;
 		  Node n1 = ngr.addNode();
 		  Node n2 = ngr.addNode();
 		  Node n3 = ngr.addNode();
 		  Node n4 = ngr.addNode();
 		  Node n5 = ngr.addNode();
 		  Node n6 = ngr.addNode();
 		  Node n7 = ngr.addNode();
 		  Arc a1 = ngr.addArc(n1, n2);
 		  Arc a2 = ngr.addArc(n1, n3);
 		  Arc a3 = ngr.addArc(n2, n4);
 		  Arc a4 = ngr.addArc(n3, n4);
 		  Arc a5 = ngr.addArc(n3, n2);
 		  Arc a6 = ngr.addArc(n5, n3);
 		  Arc a7 = ngr.addArc(n5, n6);
 		  Arc a8 = ngr.addArc(n6, n7);
 		  Arc a9 = ngr.addArc(n7, n5);
 		  Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
 		  ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
 		  Digraph::NodeMap<int> ns(ngr, 0);
 		  nl2[a7] =  1000;
 		  nl2[a8] = -1000;
 		  ns[n1] =  100;
 		  ns[n4] = -100;
 		  nc[a1] =  100;
 		  nc[a2] =   30;
 		  nc[a3] =   20;
 		  nc[a4] =   80;
 		  nc[a5] =   50;
 		  nc[a6] =   10;
 		  nc[a7] =   80;
 		  nc[a8] =   30;
 		  nc[a9] = -120;
 		  // A. Test NetworkSimplex with the default pivot rule
+		  {
 		    NetworkSimplex<Digraph> mcf(gr);
 		    // Check the equality form
 		    mcf.upperMap(u).costMap(c);
 		    checkMcf(mcf, mcf.supplyMap(s1).run(),
 		             gr, l1, u, c, s1, true,  5240, "#A1");
+		             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
 		    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
 		             gr, l1, u, c, s2, true,  7620, "#A2");
+		             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
 		    mcf.lowerMap(l2);
 		    checkMcf(mcf, mcf.supplyMap(s1).run(),
 		             gr, l2, u, c, s1, true,  5970, "#A3");
+		             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
 		    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
 		             gr, l2, u, c, s2, true,  8010, "#A4");
+		             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
 		    mcf.reset();
 		    checkMcf(mcf, mcf.supplyMap(s1).run(),
 		             gr, l1, cu, cc, s1, true,  74, "#A5");
+		             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
 		    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
 		             gr, l2, cu, cc, s2, true,  94, "#A6");
+		             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
 		    mcf.reset();
 		    checkMcf(mcf, mcf.run(),
 		             gr, l1, cu, cc, s3, true,   0, "#A7");
 		    checkMcf(mcf, mcf.boundMaps(l2, u).run(),
-		             gr, l2, u, cc, s3, false,   0, "#A8");
+		             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
 		    checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
 		             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
 		    mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
 		    checkMcf(mcf, mcf.run(),
 		             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
 		    // Check the GEQ form
-		    mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
+		    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
 		    checkMcf(mcf, mcf.run(),
 		             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
 		    mcf.problemType(mcf.GEQ);
 		             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
 		    mcf.supplyType(mcf.GEQ);
 		    checkMcf(mcf, mcf.lowerMap(l2).run(),
 		             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
 		    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
 		             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
 		    mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
 		    checkMcf(mcf, mcf.run(),
-		             gr, l2, u, c, s5, false,    0, "#A11", GEQ);
+		             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
 		    // Check the LEQ form
 		    mcf.reset().problemType(mcf.LEQ);
 		    mcf.upperMap(u).costMap(c).supplyMap(s5);
 		    mcf.reset().supplyType(mcf.LEQ);
 		    mcf.upperMap(u).costMap(c).supplyMap(s6);
 		    checkMcf(mcf, mcf.run(),
-		             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
+		             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
 		    checkMcf(mcf, mcf.lowerMap(l2).run(),
 		             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
 		    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
 		             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
 		    mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
 		    checkMcf(mcf, mcf.run(),
-		             gr, l2, u, c, s4, false,    0, "#A14", LEQ);
+		             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
 		    // Check negative costs
 		    NetworkSimplex<Digraph> nmcf(ngr);
 		    nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
 		    checkMcf(nmcf, nmcf.run(),
 		      ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
 		    checkMcf(nmcf, nmcf.upperMap(nu2).run(),
 		      ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
 		    nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
 		    checkMcf(nmcf, nmcf.run(),
 		      ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
+		  }
 		  // B. Test NetworkSimplex with each pivot rule
+		  {
 		    NetworkSimplex<Digraph> mcf(gr);
-		    mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
+		    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
 		    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
 		             gr, l2, u, c, s1, true,  5970, "#B1");
+		             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
 		    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
 		             gr, l2, u, c, s1, true,  5970, "#B2");
+		             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
 		    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
 		             gr, l2, u, c, s1, true,  5970, "#B3");
+		             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
 		    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
 		             gr, l2, u, c, s1, true,  5970, "#B4");
+		             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
 		    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
 		             gr, l2, u, c, s1, true,  5970, "#B5");
+		             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
+		  }
 		  return 0;
+		}

tools/dimacs-solver.cc

Show inline comments

@@ -114,18 +114,18 @@
 		      std::cerr << "GEQ supply contraints are used for NetworkSimplex\n\n";
 		    else
 		      std::cerr << "LEQ supply contraints are used for NetworkSimplex\n\n";
+		  }
 		  if (report) std::cerr << "Read the file: " << ti << '\n';
 		  ti.restart();
 		  NetworkSimplex<Digraph, Value> ns(g);
 		  ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
 		  if (sum_sup > 0) ns.problemType(ns.LEQ);
 		  ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup);
 		  if (sum_sup > 0) ns.supplyType(ns.LEQ);
 		  if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';
 		  ti.restart();
 		  bool res = ns.run();
 		  if (report) {
 		    std::cerr << "Run NetworkSimplex: " << ti << "\n\n";
 		    std::cerr << "Feasible flow: " << (res ? "found" : "not found") << '\n';
 		    if (res) std::cerr << "Min flow cost: " << ns.totalCost() << '\n';
+		  }

0 comments (0 inline)

RhodeCode

Login to your account