lemon: comparison lemon/network

equal deleted inserted replaced

-:09b5efc4ac5c
+:9c7a4c1327ef
 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.
 ///
 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,
 const GR &_graph;
 int _node_num;
 int _arc_num;
 // Parameters of the problem
-ValueArcMap *_plower;
+bool _have_lower;
-ValueArcMap *_pupper;
-CostArcMap *_pcost;
-ValueNodeMap *_psupply;
-bool _pstsup;
-Node _psource, _ptarget;
-Value _pstflow;
 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;
+ValueVector _supply;
-FlowVector _flow;
+ValueVector _flow;
 CostVector _pi;
 // Data for storing the spanning tree structure
 IntVector _parent;
 IntVector _pred;
 ///
 /// The constructor of the class.
 ///
 /// \param graph The digraph the algorithm runs on.
 NetworkSimplex(const GR& graph) :
-_graph(graph),
+_graph(graph), _node_id(graph), _arc_id(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),
 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
-/// Destructor.
+_node_num = countNodes(_graph);
-~NetworkSimplex() {
+_arc_num = countArcs(_graph);
-if (_local_flow) delete _flow_map;
+int all_node_num = _node_num + 1;
-if (_local_potential) delete _potential_map;
+int all_arc_num = _arc_num + _node_num;
+_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.
 /// of the algorithm.
 ///
 /// \return <tt>(*this)</tt>
 template <typename LowerMap>
 NetworkSimplex& lowerMap(const LowerMap& map) {
-delete _plower;
+_have_lower = true;
-_plower = new ValueArcMap(_graph);
 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.
 /// 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.
 /// 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.
 /// 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.
 /// \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;
+for (int i = 0; i != _node_num; ++i) {
-_psupply = NULL;
+_supply[i] = 0;
-_pstsup = true;
+}
-_psource = s;
+_supply[_node_id[s]] =  k;
-_ptarget = t;
+_supply[_node_id[t]] = -k;
-_pstflow = k;
 return *this;
 }
 /// \brief Set the type of the supply constraints.
 ///
 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();
 ///
 /// 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,
 /// \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 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.
+/// 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);
 ///
 ///     .supplyMap(sup).run();
 /// \endcode
 ///
 /// \return <tt>(*this)</tt>
 NetworkSimplex& reset() {
-delete _plower;
+for (int i = 0; i != _node_num; ++i) {
-delete _pupper;
+_supply[i] = 0;
-delete _pcost;
+}
-delete _psupply;
+for (int i = 0; i != _arc_num; ++i) {
-_plower = NULL;
+_lower[i] = 0;
-_pupper = NULL;
+_upper[i] = INF;
-_pcost = NULL;
+_cost[i] = 1;
-_psupply = NULL;
+}
-_pstsup = false;
+_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;
 }
 /// @}
 /// 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>
+template <typename Number>
-Value totalCost() const {
+Number totalCost() const {
-Value c = 0;
+Number c = 0;
-if (_pcost) {
+for (ArcIt a(_graph); a != INVALID; ++a) {
-for (ArcIt e(_graph); e != INVALID; ++e)
+int i = _arc_id[a];
-c += (*_flow_map)[e] * (*_pcost)[e];
+c += Number(_flow[i]) * Number(_cost[i]);
-} else {
-for (ArcIt e(_graph); e != INVALID; ++e)
-c += (*_flow_map)[e];
 }
 return c;
 }
 #ifndef DOXYGEN
 ///
 /// 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
+/// This function copies the flow value on each arc into the given
-/// the found flow.
+/// 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 {
+template <typename FlowMap>
-return *_flow_map;
+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
+/// \brief Return the potential map (the dual solution).
-/// (the dual solution).
+///
-///
+/// This function copies the potential (dual value) of each node
-/// This function returns a const reference to a node map storing
+/// into the given map.
-/// the found potentials, which form the dual solution of the
+/// The \c Cost type of the algorithm must be convertible to the
-/// \ref min_cost_flow "minimum cost flow problem".
+/// \c Value type of the map.
 ///
 /// \pre \ref run() must be called before using this function.
-const PotentialMap& potentialMap() const {
+template <typename PotentialMap>
-return *_potential_map;
+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);
+// Check the sum of supply values
-_source.resize(all_arc_num);
-_target.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);
-// Initialize node related data
-bool valid_supply = true;
 _sum_supply = 0;
-if (!_pstsup && !_psupply) {
+for (int i = 0; i != _node_num; ++i) {
-_pstsup = true;
+_sum_supply += _supply[i];
-_psource = _ptarget = NodeIt(_graph);
+}
-_pstflow = 0;
+if ( !(_stype == GEQ && _sum_supply <= 0 ||
-}
+_stype == LEQ && _sum_supply >= 0) ) return false;
-if (_psupply) {
-int i = 0;
+// Remove non-zero lower bounds
-for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+if (_have_lower) {
-_node_id[n] = i;
+for (int i = 0; i != _arc_num; ++i) {
-_supply[i] = (*_psupply)[n];
+Value c = _lower[i];
-_sum_supply += _supply[i];
+if (c >= 0) {
-}
+_cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
-valid_supply = (_stype == GEQ && _sum_supply <= 0) ||
+} else {
-(_stype == LEQ && _sum_supply >= 0);
+_cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
+}
+_supply[_source[i]] -= c;
+_supply[_target[i]] += c;
+}
 } else {
-int i = 0;
+for (int i = 0; i != _arc_num; ++i) {
-for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+_cap[i] = _upper[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;
 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] = _sum_supply < 0 ? -ART_COST : ART_COST;
-_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;
-}
-}
-}
 // 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];
 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 if (_flow[e] != 0) return INFEASIBLE;
 }
 }
-// Copy flow values to _flow_map
+// Transform the solution and the supply map to the original form
-if (_plower) {
+if (_have_lower) {
 for (int i = 0; i != _arc_num; ++i) {
-Arc e = _arc_ref[i];
+Value c = _lower[i];
-_flow_map->set(e, (*_plower)[e] + _flow[i]);
+if (c != 0) {
-}
+_flow[i] += c;
-} else {
+_supply[_source[i]] += c;
-for (int i = 0; i != _arc_num; ++i) {
+_supply[_target[i]] -= c;
-_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;
 }

changeset 689	111698359429
parent 688	756a5ec551c8
child 690	f3792d5bb294