LEMON/LEMON-official Changeset - r689:111698359429

Repositories » LEMON » LEMON-official

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

Changeset - r689:111698359429

« 688:756a5e

690:f3792d »

r689:111698359429

Wed, 29 Apr 2009 16:54:27

[Not Reviewed]

0 comments (0 inline)

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

Less map copying in NetworkSimplex (#234) - The graph is copied in the constructor instead of the init() function. It must not be modified after the class is constructed. - The maps are copied once (instead of twice). - Remove FlowMap, PotentialMap typedefs and flowMap(), pontentialMap() setter functions. - flowMap() and potentialMap() query functions copy the values into the given map (reference) instead of returning a const reference to a previously constructed map.

0 2 0

default

2 files changed with 184 insertions and 318 deletions:

lemon/network_simplex.h

155

286

test/min_cost_flow_test.cc

↑ Collapse diff ↑

lemon/network_simplex.h

Show inline comments

@@ -69,27 +69,16 @@
 		  /// by default. For more information see \ref PivotRule.
 		  template <typename GR, typename V = int, typename C = V>
 		  class NetworkSimplex
+		  {
 		  public:
-		    /// The flow type of the algorithm
+		    /// The type of the flow amounts, capacity bounds and supply values
 		    typedef V Value;
-		    /// The cost type of the algorithm
+		    /// The type of the arc costs
 		    typedef C Cost;
 		#ifdef DOXYGEN
 		    /// The type of the flow map
 		    typedef GR::ArcMap<Value> FlowMap;
 		    /// The type of the potential map
 		    typedef GR::NodeMap<Cost> PotentialMap;
 		#else
 		    /// The type of the flow map
 		    typedef typename GR::template ArcMap<Value> FlowMap;
 		    /// The type of the potential map
 		    typedef typename GR::template NodeMap<Cost> PotentialMap;
 		#endif
 		  public:
 		    /// \brief Problem type constants for the \c run() function.
 		    ///
 		    /// Enum type containing the problem type constants that can be
@@ -203,21 +192,17 @@
 		    };
 		  private:
 		    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 		    typedef typename GR::template ArcMap<Value> ValueArcMap;
 		    typedef typename GR::template ArcMap<Cost> CostArcMap;
 		    typedef typename GR::template NodeMap<Value> ValueNodeMap;
 		    typedef std::vector<Arc> ArcVector;
 		    typedef std::vector<Node> NodeVector;
 		    typedef std::vector<int> IntVector;
 		    typedef std::vector<bool> BoolVector;
-		    typedef std::vector<Value> FlowVector;
+		    typedef std::vector<Value> ValueVector;
 		    typedef std::vector<Cost> CostVector;
 		    // State constants for arcs
 		    enum ArcStateEnum {
 		      STATE_UPPER = -1,
 		      STATE_TREE  =  0,
@@ -229,40 +214,29 @@
 		    // Data related to the underlying digraph
 		    const GR &_graph;
 		    int _node_num;
 		    int _arc_num;
 		    // Parameters of the problem
 		    ValueArcMap *_plower;
 		    ValueArcMap *_pupper;
 		    CostArcMap *_pcost;
 		    ValueNodeMap *_psupply;
 		    bool _pstsup;
 		    Node _psource, _ptarget;
-		    Value _pstflow;
+		    bool _have_lower;
 		    SupplyType _stype;
 		    Value _sum_supply;
 		    // Result maps
 		    FlowMap *_flow_map;
 		    PotentialMap *_potential_map;
 		    bool _local_flow;
 		    bool _local_potential;
 		    // Data structures for storing the digraph
 		    IntNodeMap _node_id;
-		    ArcVector _arc_ref;
+		    IntArcMap _arc_id;
 		    IntVector _source;
 		    IntVector _target;
 		    // Node and arc data
-		    FlowVector _cap;
+		    ValueVector _lower;
 		    ValueVector _upper;
 		    ValueVector _cap;
 		    CostVector _cost;
 		    FlowVector _supply;
 		    FlowVector _flow;
 		    ValueVector _supply;
 		    ValueVector _flow;
 		    CostVector _pi;
 		    // Data for storing the spanning tree structure
 		    IntVector _parent;
 		    IntVector _pred;
 		    IntVector _thread;
@@ -686,33 +660,74 @@
 		    /// \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), _stype(GEQ),
 		      _flow_map(NULL), _potential_map(NULL),
 		      _local_flow(false), _local_potential(false),
 		      _node_id(graph),
 		      _graph(graph), _node_id(graph), _arc_id(graph),
 		      INF(std::numeric_limits<Value>::has_infinity ?
 		          std::numeric_limits<Value>::infinity() :
 		          std::numeric_limits<Value>::max())
+		    {
 		      // Check the value types
 		      LEMON_ASSERT(std::numeric_limits<Value>::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");
+		    }
 		      // Resize vectors
 		      _node_num = countNodes(_graph);
 		      _arc_num = countArcs(_graph);
 		      int all_node_num = _node_num + 1;
 		      int all_arc_num = _arc_num + _node_num;
 		    /// Destructor.
 		    ~NetworkSimplex() {
 		      if (_local_flow) delete _flow_map;
 		      if (_local_potential) delete _potential_map;
 		      _source.resize(all_arc_num);
 		      _target.resize(all_arc_num);
 		      _lower.resize(all_arc_num);
 		      _upper.resize(all_arc_num);
 		      _cap.resize(all_arc_num);
 		      _cost.resize(all_arc_num);
 		      _supply.resize(all_node_num);
 		      _flow.resize(all_arc_num);
 		      _pi.resize(all_node_num);
 		      _parent.resize(all_node_num);
 		      _pred.resize(all_node_num);
 		      _forward.resize(all_node_num);
 		      _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);
 		      // Copy the graph (store the arcs in a mixed order)
 		      int i = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
 		        _node_id[n] = i;
+		      }
 		      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
 		      i = 0;
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        _arc_id[a] = i;
 		        _source[i] = _node_id[_graph.source(a)];
 		        _target[i] = _node_id[_graph.target(a)];
 		        if ((i += k) >= _arc_num) i = (i % k) + 1;
+		      }
 		      // Initialize maps
 		      for (int i = 0; i != _node_num; ++i) {
 		        _supply[i] = 0;
+		      }
 		      for (int i = 0; i != _arc_num; ++i) {
 		        _lower[i] = 0;
 		        _upper[i] = INF;
 		        _cost[i] = 1;
+		      }
 		      _have_lower = false;
 		      _stype = GEQ;
+		    }
 		    /// \name Parameters
 		    /// The parameters of the algorithm can be specified using these
 		    /// functions.
@@ -728,16 +743,15 @@
 		    /// Its \c Value type must be convertible to the \c Value type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template <typename LowerMap>
 		    NetworkSimplex& lowerMap(const LowerMap& map) {
 		      delete _plower;
 		      _plower = new ValueArcMap(_graph);
 		      _have_lower = true;
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
-		        (*_plower)[a] = map[a];
+		        _lower[_arc_id[a]] = map[a];
+		      }
 		      return *this;
+		    }
 		    /// \brief Set the upper bounds (capacities) on the arcs.
 		    ///
@@ -750,16 +764,14 @@
 		    /// Its \c Value type must be convertible to the \c Value type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template<typename UpperMap>
 		    NetworkSimplex& upperMap(const UpperMap& map) {
 		      delete _pupper;
 		      _pupper = new ValueArcMap(_graph);
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
-		        (*_pupper)[a] = map[a];
+		        _upper[_arc_id[a]] = map[a];
+		      }
 		      return *this;
+		    }
 		    /// \brief Set the costs of the arcs.
 		    ///
@@ -771,16 +783,14 @@
 		    /// Its \c Value type must be convertible to the \c Cost type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    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];
+		        _cost[_arc_id[a]] = map[a];
+		      }
 		      return *this;
+		    }
 		    /// \brief Set the supply values of the nodes.
 		    ///
@@ -793,17 +803,14 @@
 		    /// Its \c Value type must be convertible to the \c Value type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template<typename SupplyMap>
 		    NetworkSimplex& supplyMap(const SupplyMap& map) {
 		      delete _psupply;
 		      _pstsup = false;
 		      _psupply = new ValueNodeMap(_graph);
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
-		        (*_psupply)[n] = map[n];
+		        _supply[_node_id[n]] = map[n];
+		      }
 		      return *this;
+		    }
 		    /// \brief Set single source and target nodes and a supply value.
 		    ///
@@ -821,18 +828,17 @@
 		    /// \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, Value k) {
 		      delete _psupply;
 		      _psupply = NULL;
 		      _pstsup = true;
 		      _psource = s;
 		      _ptarget = t;
 		      _pstflow = k;
 		      for (int i = 0; i != _node_num; ++i) {
 		        _supply[i] = 0;
+		      }
 		      _supply[_node_id[s]] =  k;
 		      _supply[_node_id[t]] = -k;
 		      return *this;
+		    }
 		    /// \brief Set the type of the supply constraints.
 		    ///
 		    /// This function sets the type of the supply/demand constraints.
@@ -844,71 +850,38 @@
 		    /// \return <tt>(*this)</tt>
 		    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
 		    /// automatically allocated map, of course.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    NetworkSimplex& flowMap(FlowMap& map) {
 		      if (_local_flow) {
 		        delete _flow_map;
 		        _local_flow = false;
+		      }
 		      _flow_map = &map;
 		      return *this;
+		    }
 		    /// \brief Set the potential map.
 		    ///
 		    /// This function sets the potential map, which is used for storing
 		    /// the dual solution.
 		    /// If it is not used before calling \ref run(), an instance will
 		    /// be allocated automatically. The destructor deallocates this
 		    /// automatically allocated map, of course.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    NetworkSimplex& potentialMap(PotentialMap& map) {
 		      if (_local_potential) {
 		        delete _potential_map;
 		        _local_potential = false;
+		      }
 		      _potential_map = &map;
 		      return *this;
+		    }
 		    /// @}
 		    /// \name Execution Control
 		    /// 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 costMap(), \ref supplyMap(), \ref stSupply(),
-		    /// \ref supplyType(), \ref flowMap() and \ref potentialMap().
 		    /// \ref supplyType().
 		    /// For example,
 		    /// \code
 		    ///   NetworkSimplex<ListDigraph> ns(graph);
 		    ///   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.
 		    /// However the underlying digraph must not be modified after this
 		    /// class have been constructed, since it copies and extends the graph.
 		    ///
 		    /// \param pivot_rule The pivot rule that will be used during the
 		    /// algorithm. For more information see \ref PivotRule.
 		    ///
 		    /// \return \c INFEASIBLE if no feasible flow exists,
 		    /// \n \c OPTIMAL if the problem has optimal solution
@@ -925,18 +898,19 @@
+		    }
 		    /// \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 costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(),
 		    /// \ref flowMap() and \ref potentialMap().
 		    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
 		    ///
 		    /// 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.
 		    /// However the underlying digraph must not be modified after this
 		    /// class have been constructed, since it copies and extends the graph.
 		    ///
 		    /// For example,
 		    /// \code
 		    ///   NetworkSimplex<ListDigraph> ns(graph);
 		    ///
 		    ///   // First run
@@ -954,29 +928,22 @@
 		    ///   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;
+		      for (int i = 0; i != _node_num; ++i) {
 		        _supply[i] = 0;
+		      }
 		      for (int i = 0; i != _arc_num; ++i) {
 		        _lower[i] = 0;
 		        _upper[i] = INF;
 		        _cost[i] = 1;
+		      }
 		      _have_lower = false;
 		      _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;
+		    }
 		    /// @}
 		    /// \name Query Functions
@@ -998,21 +965,18 @@
 		    /// \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 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];
+		    template <typename Number>
 		    Number totalCost() const {
 		      Number c = 0;
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        int i = _arc_id[a];
 		        c += Number(_flow[i]) * Number(_cost[i]);
+		      }
 		      return c;
+		    }
 		#ifndef DOXYGEN
 		    Cost totalCost() const {
@@ -1023,116 +987,87 @@
 		    /// \brief Return the flow on the given arc.
 		    ///
 		    /// This function returns the flow on the given arc.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    Value flow(const Arc& a) const {
-		      return (*_flow_map)[a];
+		      return _flow[_arc_id[a]];
+		    }
-		    /// \brief Return a const reference to the flow map.
+		    /// \brief Return the flow map (the primal solution).
 		    ///
 		    /// This function returns a const reference to an arc map storing
 		    /// the found flow.
 		    /// This function copies the flow value on each arc into the given
 		    /// map. The \c Value type of the algorithm must be convertible to
 		    /// the \c Value type of the map.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    const FlowMap& flowMap() const {
 		      return *_flow_map;
 		    template <typename FlowMap>
 		    void flowMap(FlowMap &map) const {
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        map.set(a, _flow[_arc_id[a]]);
+		      }
+		    }
 		    /// \brief Return the potential (dual value) of the given node.
 		    ///
 		    /// This function returns the potential (dual value) of the
 		    /// given node.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    Cost potential(const Node& n) const {
-		      return (*_potential_map)[n];
+		      return _pi[_node_id[n]];
+		    }
 		    /// \brief Return a const reference to the potential map
 		    /// (the dual solution).
 		    /// \brief Return 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".
+		    /// This function copies the potential (dual value) of each node
 		    /// into the given map.
 		    /// The \c Cost type of the algorithm must be convertible to the
 		    /// \c Value type of the map.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    const PotentialMap& potentialMap() const {
 		      return *_potential_map;
 		    template <typename PotentialMap>
 		    void potentialMap(PotentialMap &map) const {
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        map.set(n, _pi[_node_id[n]]);
+		      }
+		    }
 		    /// @}
 		  private:
 		    // Initialize internal data structures
 		    bool init() {
 		      // Initialize result maps
 		      if (!_flow_map) {
 		        _flow_map = new FlowMap(_graph);
 		        _local_flow = true;
+		      }
 		      if (!_potential_map) {
 		        _potential_map = new PotentialMap(_graph);
 		        _local_potential = true;
+		      }
 		      // Initialize vectors
 		      _node_num = countNodes(_graph);
 		      _arc_num = countArcs(_graph);
 		      int all_node_num = _node_num + 1;
 		      int all_arc_num = _arc_num + _node_num;
 		      if (_node_num == 0) return false;
 		      _arc_ref.resize(_arc_num);
 		      _source.resize(all_arc_num);
-		      _target.resize(all_arc_num);
+		      // Check the sum of supply values
 		      _sum_supply = 0;
 		      for (int i = 0; i != _node_num; ++i) {
 		        _sum_supply += _supply[i];
+		      }
 		      if ( !(_stype == GEQ && _sum_supply <= 0 ||
 		             _stype == LEQ && _sum_supply >= 0) ) return false;
 		      _cap.resize(all_arc_num);
 		      _cost.resize(all_arc_num);
 		      _supply.resize(all_node_num);
 		      _flow.resize(all_arc_num);
 		      _pi.resize(all_node_num);
 		      _parent.resize(all_node_num);
 		      _pred.resize(all_node_num);
 		      _forward.resize(all_node_num);
 		      _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;
 		      _sum_supply = 0;
 		      if (!_pstsup && !_psupply) {
 		        _pstsup = true;
 		        _psource = _ptarget = NodeIt(_graph);
 		        _pstflow = 0;
 		      // Remove non-zero lower bounds
 		      if (_have_lower) {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Value c = _lower[i];
 		          if (c >= 0) {
 		            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
 		          } else {
 		            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
+		          }
 		          _supply[_source[i]] -= c;
 		          _supply[_target[i]] += c;
+		        }
 		      } else {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          _cap[i] = _upper[i];
+		        }
+		      }
 		      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];
+		        }
 		        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;
 		      // Initialize artifical cost
 		      Cost ART_COST;
 		      if (std::numeric_limits<Cost>::is_exact) {
 		        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
 		      } else {
@@ -1140,99 +1075,37 @@
 		        for (int i = 0; i != _arc_num; ++i) {
 		          if (_cost[i] > ART_COST) ART_COST = _cost[i];
+		        }
 		        ART_COST = (ART_COST + 1) * _node_num;
+		      }
 		      // Initialize arc maps
 		      for (int i = 0; i != _arc_num; ++i) {
 		        _flow[i] = 0;
 		        _state[i] = STATE_LOWER;
+		      }
 		      // 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;
+		      _succ_num[_root] = _node_num + 1;
 		      _last_succ[_root] = _root - 1;
 		      _supply[_root] = -_sum_supply;
 		      if (_sum_supply < 0) {
 		        _pi[_root] = -ART_COST;
 		      } else {
 		        _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;
+		      }
 		      // Initialize arc maps
 		      if (_pupper && _pcost) {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Arc e = _arc_ref[i];
 		          _source[i] = _node_id[_graph.source(e)];
 		          _target[i] = _node_id[_graph.target(e)];
 		          _cap[i] = (*_pupper)[e];
 		          _cost[i] = (*_pcost)[e];
 		          _flow[i] = 0;
 		          _state[i] = STATE_LOWER;
+		        }
 		      } else {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Arc e = _arc_ref[i];
 		          _source[i] = _node_id[_graph.source(e)];
 		          _target[i] = _node_id[_graph.target(e)];
 		          _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;
+		        }
 		        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) {
 		          Value c = (*_plower)[_arc_ref[i]];
 		          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;
+		          }
+		        }
+		      }
+		      _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
 		      // Add artificial arcs and initialize the spanning tree data structure
 		      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 		        _parent[u] = _root;
 		        _pred[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;
 		        _state[e] = STATE_TREE;
 		        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
 		          _flow[e] = _supply[u];
 		          _forward[u] = true;
@@ -1514,26 +1387,22 @@
 		      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) {
 		      // Transform the solution and the supply map to the original form
 		      if (_have_lower) {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Arc e = _arc_ref[i];
 		          _flow_map->set(e, (*_plower)[e] + _flow[i]);
 		          Value c = _lower[i];
 		          if (c != 0) {
 		            _flow[i] += c;
 		            _supply[_source[i]] += c;
 		            _supply[_target[i]] -= c;
+		          }
+		        }
 		      } else {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          _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 OPTIMAL;
+		    }
 		  }; //class NetworkSimplex

test/min_cost_flow_test.cc

Show inline comments

@@ -81,71 +81,63 @@
 		  EQ,
 		  GEQ,
 		  LEQ
 		};
 		// Check the interface of an MCF algorithm
-		template <typename GR, typename Flow, typename Cost>
+		template <typename GR, typename Value, typename Cost>
 		class McfClassConcept
+		{
 		public:
 		  template <typename MCF>
 		  struct Constraints {
 		    void constraints() {
 		      checkConcept<concepts::Digraph, GR>();
 		      MCF mcf(g);
 		      const MCF& const_mcf = mcf;
 		      b = mcf.reset()
 		             .lowerMap(lower)
 		             .upperMap(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();
 		      c = const_mcf.totalCost();
-		      double x = const_mcf.template totalCost<double>();
 		      x = const_mcf.template totalCost<double>();
 		      v = const_mcf.flow(a);
 		      c = const_mcf.potential(n);
 		      v = const_mcf.INF;
 		      ignore_unused_variable_warning(fm);
 		      ignore_unused_variable_warning(pm);
 		      ignore_unused_variable_warning(x);
 		      const_mcf.flowMap(fm);
 		      const_mcf.potentialMap(pm);
+		    }
 		    typedef typename GR::Node Node;
 		    typedef typename GR::Arc Arc;
 		    typedef concepts::ReadMap<Node, Flow> NM;
 		    typedef concepts::ReadMap<Arc, Flow> FAM;
 		    typedef concepts::ReadMap<Node, Value> NM;
 		    typedef concepts::ReadMap<Arc, Value> VAM;
 		    typedef concepts::ReadMap<Arc, Cost> CAM;
 		    typedef concepts::WriteMap<Arc, Value> FlowMap;
 		    typedef concepts::WriteMap<Node, Cost> PotMap;
 		    const GR &g;
 		    const FAM &lower;
 		    const FAM &upper;
 		    const VAM &lower;
 		    const VAM &upper;
 		    const CAM &cost;
 		    const NM &sup;
 		    const Node &n;
 		    const Arc &a;
 		    const Flow &k;
 		    Flow v;
-		    Cost c;
+		    const Value &k;
 		    FlowMap fm;
 		    PotMap pm;
 		    bool b;
 		    typename MCF::FlowMap &flow;
-		    typename MCF::PotentialMap &pot;
+		    double x;
 		    typename MCF::Value v;
 		    typename MCF::Cost c;
 		  };
 		};
 		// Check the feasibility of the given flow (primal soluiton)
@@ -218,30 +210,35 @@
 		               PT result, bool optimal, typename CM::Value total,
 		               const std::string &test_id = "",
 		               SupplyType type = EQ )
+		{
 		  check(mcf_result == result, "Wrong result " + test_id);
 		  if (optimal) {
-		    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
+		    typename GR::template ArcMap<typename SM::Value> flow(gr);
 		    typename GR::template NodeMap<typename CM::Value> pi(gr);
 		    mcf.flowMap(flow);
 		    mcf.potentialMap(pi);
 		    check(checkFlow(gr, lower, upper, supply, flow, 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()),
 		    check(checkPotential(gr, lower, upper, cost, supply, flow, pi),
 		          "Wrong potentials " + test_id);
+		  }
+		}
 		int main()
+		{
 		  // Check the interfaces
+		  {
 		    typedef int Flow;
 		    typedef int Cost;
 		    typedef concepts::Digraph GR;
 		    checkConcept< McfClassConcept<GR, Flow, Cost>,
 		                  NetworkSimplex<GR, Flow, Cost> >();
 		    checkConcept< McfClassConcept<GR, int, int>,
 		                  NetworkSimplex<GR> >();
 		    checkConcept< McfClassConcept<GR, double, double>,
 		                  NetworkSimplex<GR, double> >();
 		    checkConcept< McfClassConcept<GR, int, double>,
 		                  NetworkSimplex<GR, int, double> >();
+		  }
 		  // Run various MCF tests
 		  typedef ListDigraph Digraph;
 		  DIGRAPH_TYPEDEFS(ListDigraph);

0 comments (0 inline)

RhodeCode

Login to your account