lemon-0.x: comparison lemon/network

equal deleted inserted replaced

-:d5c822bd99ec
+:13ab2cd16858
 #include <vector>
 #include <limits>
 #include <algorithm>
-#include <lemon/graph_adaptor.h>
 #include <lemon/graph_utils.h>
-#include <lemon/smart_graph.h>
 #include <lemon/math.h>
 namespace lemon {
 /// \addtogroup min_cost_flow
 typename CapacityMap = typename Graph::template EdgeMap<int>,
 typename CostMap = typename Graph::template EdgeMap<int>,
 typename SupplyMap = typename Graph::template NodeMap<int> >
 class NetworkSimplex
 {
+GRAPH_TYPEDEFS(typename Graph);
 typedef typename CapacityMap::Value Capacity;
 typedef typename CostMap::Value Cost;
 typedef typename SupplyMap::Value Supply;
-typedef SmartGraph SGraph;
-GRAPH_TYPEDEFS(typename SGraph);
-typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
-typedef typename SGraph::template EdgeMap<Cost> SCostMap;
-typedef typename SGraph::template NodeMap<Supply> SSupplyMap;
-typedef typename SGraph::template NodeMap<Cost> SPotentialMap;
-typedef typename SGraph::template NodeMap<int> IntNodeMap;
-typedef typename SGraph::template NodeMap<bool> BoolNodeMap;
-typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
-typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
-typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
-typedef typename SGraph::template EdgeMap<bool> BoolEdgeMap;
-typedef typename Graph::template NodeMap<Node> NodeRefMap;
-typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
 typedef std::vector<Edge> EdgeVector;
+typedef std::vector<Node> NodeVector;
+typedef std::vector<int> IntVector;
+typedef std::vector<bool> BoolVector;
+typedef std::vector<Capacity> CapacityVector;
+typedef std::vector<Cost> CostVector;
+typedef std::vector<Supply> SupplyVector;
+typedef typename Graph::template NodeMap<int> IntNodeMap;
 public:
 /// The type of the flow map.
 typedef typename Graph::template EdgeMap<Capacity> FlowMap;
 ALTERING_LIST_PIVOT
 };
 private:
-/// \brief Map adaptor class for handling reduced edge costs.
-///
-/// Map adaptor class for handling reduced edge costs.
-class ReducedCostMap : public MapBase<Edge, Cost>
-{
-private:
-const SGraph &_gr;
-const SCostMap &_cost_map;
-const SPotentialMap &_pot_map;
-public:
-///\e
-ReducedCostMap( const SGraph &gr,
-const SCostMap &cost_map,
-const SPotentialMap &pot_map ) :
-_gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
-///\e
-Cost operator[](const Edge &e) const {
-return _cost_map[e] + _pot_map[_gr.source(e)]
-- _pot_map[_gr.target(e)];
-}
-}; //class ReducedCostMap
-private:
 /// \brief Implementation of the "First Eligible" pivot rule for the
 /// \ref NetworkSimplex "network simplex" algorithm.
 ///
 /// This class implements the "First Eligible" pivot rule
 /// for the \ref NetworkSimplex "network simplex" algorithm.
 class FirstEligiblePivotRule
 {
 private:
 // References to the NetworkSimplex class
-NetworkSimplex &_ns;
+const EdgeVector &_edge;
-EdgeVector &_edges;
+const IntVector  &_source;
+const IntVector  &_target;
+const CostVector &_cost;
+const IntVector  &_state;
+const CostVector &_pi;
+int &_in_edge;
+int _edge_num;
+// Pivot rule data
 int _next_edge;
 public:
 /// Constructor
-FirstEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+FirstEligiblePivotRule(NetworkSimplex &ns) :
-_ns(ns), _edges(edges), _next_edge(0) {}
+_edge(ns._edge), _source(ns._source), _target(ns._target),
+_cost(ns._cost), _state(ns._state), _pi(ns._pi),
+_in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
+{}
 /// Find next entering edge
 bool findEnteringEdge() {
-Edge e;
+Cost c;
-for (int i = _next_edge; i < int(_edges.size()); ++i) {
+for (int e = _next_edge; e < _edge_num; ++e) {
-e = _edges[i];
+c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-if (_ns._state[e] * _ns._red_cost[e] < 0) {
+if (c < 0) {
-_ns._in_edge = e;
+_in_edge = e;
-_next_edge = i + 1;
+_next_edge = e + 1;
 return true;
 }
 }
-for (int i = 0; i < _next_edge; ++i) {
+for (int e = 0; e < _next_edge; ++e) {
-e = _edges[i];
+c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-if (_ns._state[e] * _ns._red_cost[e] < 0) {
+if (c < 0) {
-_ns._in_edge = e;
+_in_edge = e;
-_next_edge = i + 1;
+_next_edge = e + 1;
 return true;
 }
 }
 return false;
 }
 }; //class FirstEligiblePivotRule
 /// \brief Implementation of the "Best Eligible" pivot rule for the
 /// \ref NetworkSimplex "network simplex" algorithm.
 ///
 /// This class implements the "Best Eligible" pivot rule
 class BestEligiblePivotRule
 {
 private:
 // References to the NetworkSimplex class
-NetworkSimplex &_ns;
+const EdgeVector &_edge;
-EdgeVector &_edges;
+const IntVector  &_source;
+const IntVector  &_target;
+const CostVector &_cost;
+const IntVector  &_state;
+const CostVector &_pi;
+int &_in_edge;
+int _edge_num;
 public:
 /// Constructor
-BestEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+BestEligiblePivotRule(NetworkSimplex &ns) :
-_ns(ns), _edges(edges) {}
+_edge(ns._edge), _source(ns._source), _target(ns._target),
+_cost(ns._cost), _state(ns._state), _pi(ns._pi),
+_in_edge(ns._in_edge), _edge_num(ns._edge_num)
+{}
 /// Find next entering edge
 bool findEnteringEdge() {
-Cost min = 0;
+Cost c, min = 0;
-Edge e;
+for (int e = 0; e < _edge_num; ++e) {
-for (int i = 0; i < int(_edges.size()); ++i) {
+c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-e = _edges[i];
+if (c < min) {
-if (_ns._state[e] * _ns._red_cost[e] < min) {
+min = c;
-min = _ns._state[e] * _ns._red_cost[e];
+_in_edge = e;
-_ns._in_edge = e;
 }
 }
 return min < 0;
 }
 }; //class BestEligiblePivotRule
 /// \brief Implementation of the "Block Search" pivot rule for the
 /// \ref NetworkSimplex "network simplex" algorithm.
 ///
 /// This class implements the "Block Search" pivot rule
 class BlockSearchPivotRule
 {
 private:
 // References to the NetworkSimplex class
-NetworkSimplex &_ns;
+const EdgeVector &_edge;
-EdgeVector &_edges;
+const IntVector  &_source;
+const IntVector  &_target;
+const CostVector &_cost;
+const IntVector  &_state;
+const CostVector &_pi;
+int &_in_edge;
+int _edge_num;
+// Pivot rule data
 int _block_size;
-int _next_edge, _min_edge;
+int _next_edge;
 public:
 /// Constructor
-BlockSearchPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+BlockSearchPivotRule(NetworkSimplex &ns) :
-_ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
+_edge(ns._edge), _source(ns._source), _target(ns._target),
+_cost(ns._cost), _state(ns._state), _pi(ns._pi),
+_in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
 {
 // The main parameters of the pivot rule
 const double BLOCK_SIZE_FACTOR = 2.0;
 const int MIN_BLOCK_SIZE = 10;
-_block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
+_block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)),
 MIN_BLOCK_SIZE );
 }
 /// Find next entering edge
 bool findEnteringEdge() {
-Cost curr, min = 0;
+Cost c, min = 0;
-Edge e;
 int cnt = _block_size;
-int i;
+int e, min_edge = _next_edge;
-for (i = _next_edge; i < int(_edges.size()); ++i) {
+for (e = _next_edge; e < _edge_num; ++e) {
-e = _edges[i];
+c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
+if (c < min) {
-min = curr;
+min = c;
-_min_edge = i;
+min_edge = e;
 }
 if (--cnt == 0) {
 if (min < 0) break;
 cnt = _block_size;
 }
 }
 if (min == 0 || cnt > 0) {
-for (i = 0; i < _next_edge; ++i) {
+for (e = 0; e < _next_edge; ++e) {
-e = _edges[i];
+c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
+if (c < min) {
-min = curr;
+min = c;
-_min_edge = i;
+min_edge = e;
 }
 if (--cnt == 0) {
 if (min < 0) break;
 cnt = _block_size;
 }
 }
 }
 if (min >= 0) return false;
-_ns._in_edge = _edges[_min_edge];
+_in_edge = min_edge;
-_next_edge = i;
+_next_edge = e;
 return true;
 }
 }; //class BlockSearchPivotRule
 /// \brief Implementation of the "Candidate List" pivot rule for the
 /// \ref NetworkSimplex "network simplex" algorithm.
 ///
 /// This class implements the "Candidate List" pivot rule
 class CandidateListPivotRule
 {
 private:
 // References to the NetworkSimplex class
-NetworkSimplex &_ns;
+const EdgeVector &_edge;
-EdgeVector &_edges;
+const IntVector  &_source;
+const IntVector  &_target;
-EdgeVector _candidates;
+const CostVector &_cost;
+const IntVector  &_state;
+const CostVector &_pi;
+int &_in_edge;
+int _edge_num;
+// Pivot rule data
+IntVector _candidates;
 int _list_length, _minor_limit;
 int _curr_length, _minor_count;
-int _next_edge, _min_edge;
+int _next_edge;
 public:
 /// Constructor
-CandidateListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+CandidateListPivotRule(NetworkSimplex &ns) :
-_ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
+_edge(ns._edge), _source(ns._source), _target(ns._target),
+_cost(ns._cost), _state(ns._state), _pi(ns._pi),
+_in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
 {
 // The main parameters of the pivot rule
 const double LIST_LENGTH_FACTOR = 1.0;
 const int MIN_LIST_LENGTH = 10;
 const double MINOR_LIMIT_FACTOR = 0.1;
 const int MIN_MINOR_LIMIT = 3;
-_list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edges.size())),
+_list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edge_num)),
 MIN_LIST_LENGTH );
 _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
 MIN_MINOR_LIMIT );
 _curr_length = _minor_count = 0;
 _candidates.resize(_list_length);
 }
 /// Find next entering edge
 bool findEnteringEdge() {
-Cost min, curr;
+Cost min, c;
+int e, min_edge = _next_edge;
 if (_curr_length > 0 && _minor_count < _minor_limit) {
 // Minor iteration: select the best eligible edge from the
 // current candidate list
 ++_minor_count;
-Edge e;
 min = 0;
 for (int i = 0; i < _curr_length; ++i) {
 e = _candidates[i];
-curr = _ns._state[e] * _ns._red_cost[e];
+c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-if (curr < min) {
+if (c < min) {
-min = curr;
+min = c;
-_ns._in_edge = e;
+min_edge = e;
 }
-if (curr >= 0) {
+if (c >= 0) {
 _candidates[i--] = _candidates[--_curr_length];
 }
 }
-if (min < 0) return true;
+if (min < 0) {
+_in_edge = min_edge;
+return true;
+}
 }
 // Major iteration: build a new candidate list
-Edge e;
 min = 0;
 _curr_length = 0;
-int i;
+for (e = _next_edge; e < _edge_num; ++e) {
-for (i = _next_edge; i < int(_edges.size()); ++i) {
+c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-e = _edges[i];
+if (c < 0) {
-if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
 _candidates[_curr_length++] = e;
-if (curr < min) {
+if (c < min) {
-min = curr;
+min = c;
-_min_edge = i;
+min_edge = e;
 }
 if (_curr_length == _list_length) break;
 }
 }
 if (_curr_length < _list_length) {
-for (i = 0; i < _next_edge; ++i) {
+for (e = 0; e < _next_edge; ++e) {
-e = _edges[i];
+c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
+if (c < 0) {
 _candidates[_curr_length++] = e;
-if (curr < min) {
+if (c < min) {
-min = curr;
+min = c;
-_min_edge = i;
+min_edge = e;
 }
 if (_curr_length == _list_length) break;
 }
 }
 }
 if (_curr_length == 0) return false;
 _minor_count = 1;
-_ns._in_edge = _edges[_min_edge];
+_in_edge = min_edge;
-_next_edge = i;
+_next_edge = e;
 return true;
 }
 }; //class CandidateListPivotRule
 /// \brief Implementation of the "Altering Candidate List" pivot rule
 /// for the \ref NetworkSimplex "network simplex" algorithm.
 ///
 /// This class implements the "Altering Candidate List" pivot rule
 class AlteringListPivotRule
 {
 private:
 // References to the NetworkSimplex class
-NetworkSimplex &_ns;
+const EdgeVector &_edge;
-EdgeVector &_edges;
+const IntVector  &_source;
+const IntVector  &_target;
-EdgeVector _candidates;
+const CostVector &_cost;
-SCostMap _cand_cost;
+const IntVector  &_state;
+const CostVector &_pi;
+int &_in_edge;
+int _edge_num;
 int _block_size, _head_length, _curr_length;
 int _next_edge;
+IntVector _candidates;
+CostVector _cand_cost;
 // Functor class to compare edges during sort of the candidate list
 class SortFunc
 {
 private:
-const SCostMap &_map;
+const CostVector &_map;
 public:
-SortFunc(const SCostMap &map) : _map(map) {}
+SortFunc(const CostVector &map) : _map(map) {}
-bool operator()(const Edge &e1, const Edge &e2) {
+bool operator()(int left, int right) {
-return _map[e1] > _map[e2];
+return _map[left] > _map[right];
 }
 };
 SortFunc _sort_func;
 public:
 /// Constructor
-AlteringListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+AlteringListPivotRule(NetworkSimplex &ns) :
-_ns(ns), _edges(edges), _cand_cost(_ns._graph),
+_edge(ns._edge), _source(ns._source), _target(ns._target),
-_next_edge(0), _sort_func(_cand_cost)
+_cost(ns._cost), _state(ns._state), _pi(ns._pi),
+_in_edge(ns._in_edge), _edge_num(ns._edge_num),
+_next_edge(0), _cand_cost(ns._edge_num),_sort_func(_cand_cost)
 {
 // The main parameters of the pivot rule
 const double BLOCK_SIZE_FACTOR = 1.5;
 const int MIN_BLOCK_SIZE = 10;
 const double HEAD_LENGTH_FACTOR = 0.1;
 const int MIN_HEAD_LENGTH = 3;
-_block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
+_block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)),
 MIN_BLOCK_SIZE );
 _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
 MIN_HEAD_LENGTH );
 _candidates.resize(_head_length + _block_size);
 _curr_length = 0;
 }
 /// Find next entering edge
 bool findEnteringEdge() {
 // Check the current candidate list
-Edge e;
+int e;
-for (int ix = 0; ix < _curr_length; ++ix) {
+for (int i = 0; i < _curr_length; ++i) {
-e = _candidates[ix];
+e = _candidates[i];
-if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) {
+_cand_cost[e] = _state[e] *
-_candidates[ix--] = _candidates[--_curr_length];
+(_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+if (_cand_cost[e] >= 0) {
+_candidates[i--] = _candidates[--_curr_length];
 }
 }
 // Extend the list
 int cnt = _block_size;
 int last_edge = 0;
 int limit = _head_length;
-for (int i = _next_edge; i < int(_edges.size()); ++i) {
-e = _edges[i];
+for (int e = _next_edge; e < _edge_num; ++e) {
-if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
+_cand_cost[e] = _state[e] *
+(_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+if (_cand_cost[e] < 0) {
 _candidates[_curr_length++] = e;
-last_edge = i;
+last_edge = e;
 }
 if (--cnt == 0) {
 if (_curr_length > limit) break;
 limit = 0;
 cnt = _block_size;
 }
 }
 if (_curr_length <= limit) {
-for (int i = 0; i < _next_edge; ++i) {
+for (int e = 0; e < _next_edge; ++e) {
-e = _edges[i];
+_cand_cost[e] = _state[e] *
-if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
+(_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+if (_cand_cost[e] < 0) {
 _candidates[_curr_length++] = e;
-last_edge = i;
+last_edge = e;
 }
 if (--cnt == 0) {
 if (_curr_length > limit) break;
 limit = 0;
 cnt = _block_size;
 // Make heap of the candidate list (approximating a partial sort)
 make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
 _sort_func );
 // Pop the first element of the heap
-_ns._in_edge = _candidates[0];
+_in_edge = _candidates[0];
 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
 _sort_func );
 _curr_length = std::min(_head_length, _curr_length - 1);
 return true;
 }
 }; //class AlteringListPivotRule
 private:
 // State constants for edges
 STATE_LOWER =  1
 };
 private:
-// The directed graph the algorithm runs on
-SGraph _graph;
 // The original graph
-const Graph &_graph_ref;
+const Graph &_orig_graph;
 // The original lower bound map
-const LowerMap *_lower;
+const LowerMap *_orig_lower;
-// The capacity map
+// The original capacity map
-SCapacityMap _capacity;
+const CapacityMap &_orig_cap;
-// The cost map
+// The original cost map
-SCostMap _cost;
+const CostMap &_orig_cost;
-// The supply map
+// The original supply map
-SSupplyMap _supply;
+const SupplyMap *_orig_supply;
-bool _valid_supply;
+// The source node (if no supply map was given)
+Node _orig_source;
-// Edge map of the current flow
+// The target node (if no supply map was given)
-SCapacityMap _flow;
+Node _orig_target;
-// Node map of the current potentials
+// The flow value (if no supply map was given)
-SPotentialMap _potential;
+Capacity _orig_flow_value;
-// The depth node map of the spanning tree structure
+// The flow result map
-IntNodeMap _depth;
-// The parent node map of the spanning tree structure
-NodeNodeMap _parent;
-// The pred_edge node map of the spanning tree structure
-EdgeNodeMap _pred_edge;
-// The thread node map of the spanning tree structure
-NodeNodeMap _thread;
-// The forward node map of the spanning tree structure
-BoolNodeMap _forward;
-// The state edge map
-IntEdgeMap _state;
-// The artificial root node of the spanning tree structure
-Node _root;
-// The reduced cost map
-ReducedCostMap _red_cost;
-// The non-artifical edges
-EdgeVector _edges;
-// Members for handling the original graph
 FlowMap *_flow_result;
+// The potential result map
 PotentialMap *_potential_result;
+// Indicate if the flow result map is local
 bool _local_flow;
+// Indicate if the potential result map is local
 bool _local_potential;
-NodeRefMap _node_ref;
-EdgeRefMap _edge_ref;
+// The edge references
+EdgeVector _edge;
+// The node references
+NodeVector _node;
+// The node ids
+IntNodeMap _node_id;
+// The source nodes of the edges
+IntVector _source;
+// The target nodess of the edges
+IntVector _target;
+// The (modified) capacity vector
+CapacityVector _cap;
+// The cost vector
+CostVector _cost;
+// The (modified) supply vector
+CostVector _supply;
+// The current flow vector
+CapacityVector _flow;
+// The current potential vector
+CostVector _pi;
+// The number of nodes in the original graph
+int _node_num;
+// The number of edges in the original graph
+int _edge_num;
+// The depth vector of the spanning tree structure
+IntVector _depth;
+// The parent vector of the spanning tree structure
+IntVector _parent;
+// The pred_edge vector of the spanning tree structure
+IntVector _pred;
+// The thread vector of the spanning tree structure
+IntVector _thread;
+// The forward vector of the spanning tree structure
+BoolVector _forward;
+// The state vector
+IntVector _state;
+// The root node
+int _root;
 // The entering edge of the current pivot iteration
-Edge _in_edge;
+int _in_edge;
 // Temporary nodes used in the current pivot iteration
-Node join, u_in, v_in, u_out, v_out;
+int join, u_in, v_in, u_out, v_out;
-Node right, first, second, last;
+int right, first, second, last;
-Node stem, par_stem, new_stem;
+int stem, par_stem, new_stem;
 // The maximum augment amount along the found cycle in the current
 // pivot iteration
 Capacity delta;
-public :
+public:
 /// \brief General constructor (with lower bounds).
 ///
 /// General constructor (with lower bounds).
 ///
 NetworkSimplex( const Graph &graph,
 const LowerMap &lower,
 const CapacityMap &capacity,
 const CostMap &cost,
 const SupplyMap &supply ) :
-_graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
+_orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
-_cost(_graph), _supply(_graph), _flow(_graph),
+_orig_cost(cost), _orig_supply(&supply),
-_potential(_graph), _depth(_graph), _parent(_graph),
-_pred_edge(_graph), _thread(_graph), _forward(_graph),
-_state(_graph), _red_cost(_graph, _cost, _potential),
 _flow_result(NULL), _potential_result(NULL),
 _local_flow(false), _local_potential(false),
-_node_ref(graph), _edge_ref(graph)
+_node_id(graph)
-{
+{}
-// Check the sum of supply values
-Supply sum = 0;
-for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
-sum += supply[n];
-if (!(_valid_supply = sum == 0)) return;
-// Copy _graph_ref to _graph
-_graph.reserveNode(countNodes(_graph_ref) + 1);
-_graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
-copyGraph(_graph, _graph_ref)
-.edgeMap(_capacity, capacity)
-.edgeMap(_cost, cost)
-.nodeMap(_supply, supply)
-.nodeRef(_node_ref)
-.edgeRef(_edge_ref)
-.run();
-// Remove non-zero lower bounds
-for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
-if (lower[e] != 0) {
-_capacity[_edge_ref[e]] = capacity[e] - lower[e];
-_supply[_node_ref[_graph_ref.source(e)]] -= lower[e];
-_supply[_node_ref[_graph_ref.target(e)]] += lower[e];
-}
-}
-}
 /// \brief General constructor (without lower bounds).
 ///
 /// General constructor (without lower bounds).
 ///
 /// \param supply The supply values of the nodes (signed).
 NetworkSimplex( const Graph &graph,
 const CapacityMap &capacity,
 const CostMap &cost,
 const SupplyMap &supply ) :
-_graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
+_orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
-_cost(_graph), _supply(_graph), _flow(_graph),
+_orig_cost(cost), _orig_supply(&supply),
-_potential(_graph), _depth(_graph), _parent(_graph),
-_pred_edge(_graph), _thread(_graph), _forward(_graph),
-_state(_graph), _red_cost(_graph, _cost, _potential),
 _flow_result(NULL), _potential_result(NULL),
 _local_flow(false), _local_potential(false),
-_node_ref(graph), _edge_ref(graph)
+_node_id(graph)
-{
+{}
-// Check the sum of supply values
-Supply sum = 0;
-for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
-sum += supply[n];
-if (!(_valid_supply = sum == 0)) return;
-// Copy _graph_ref to _graph
-_graph.reserveNode(countNodes(_graph_ref) + 1);
-_graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
-copyGraph(_graph, _graph_ref)
-.edgeMap(_capacity, capacity)
-.edgeMap(_cost, cost)
-.nodeMap(_supply, supply)
-.nodeRef(_node_ref)
-.edgeRef(_edge_ref)
-.run();
-}
 /// \brief Simple constructor (with lower bounds).
 ///
 /// Simple constructor (with lower bounds).
 ///
 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
 NetworkSimplex( const Graph &graph,
 const LowerMap &lower,
 const CapacityMap &capacity,
 const CostMap &cost,
-typename Graph::Node s,
+Node s, Node t,
-typename Graph::Node t,
+Capacity flow_value ) :
-typename SupplyMap::Value flow_value ) :
+_orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
-_graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
+_orig_cost(cost), _orig_supply(NULL),
-_cost(_graph), _supply(_graph, 0), _flow(_graph),
+_orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
-_potential(_graph), _depth(_graph), _parent(_graph),
-_pred_edge(_graph), _thread(_graph), _forward(_graph),
-_state(_graph), _red_cost(_graph, _cost, _potential),
 _flow_result(NULL), _potential_result(NULL),
 _local_flow(false), _local_potential(false),
-_node_ref(graph), _edge_ref(graph)
+_node_id(graph)
-{
+{}
-// Copy _graph_ref to graph
-_graph.reserveNode(countNodes(_graph_ref) + 1);
-_graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
-copyGraph(_graph, _graph_ref)
-.edgeMap(_capacity, capacity)
-.edgeMap(_cost, cost)
-.nodeRef(_node_ref)
-.edgeRef(_edge_ref)
-.run();
-// Remove non-zero lower bounds
-_supply[_node_ref[s]] =  flow_value;
-_supply[_node_ref[t]] = -flow_value;
-for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
-if (lower[e] != 0) {
-_capacity[_edge_ref[e]] = capacity[e] - lower[e];
-_supply[_node_ref[_graph_ref.source(e)]] -= lower[e];
-_supply[_node_ref[_graph_ref.target(e)]] += lower[e];
-}
-}
-_valid_supply = true;
-}
 /// \brief Simple constructor (without lower bounds).
 ///
 /// Simple constructor (without lower bounds).
 ///
 /// \param flow_value 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).
 NetworkSimplex( const Graph &graph,
 const CapacityMap &capacity,
 const CostMap &cost,
-typename Graph::Node s,
+Node s, Node t,
-typename Graph::Node t,
+Capacity flow_value ) :
-typename SupplyMap::Value flow_value ) :
+_orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
-_graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
+_orig_cost(cost), _orig_supply(NULL),
-_cost(_graph), _supply(_graph, 0), _flow(_graph),
+_orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
-_potential(_graph), _depth(_graph), _parent(_graph),
-_pred_edge(_graph), _thread(_graph), _forward(_graph),
-_state(_graph), _red_cost(_graph, _cost, _potential),
 _flow_result(NULL), _potential_result(NULL),
 _local_flow(false), _local_potential(false),
-_node_ref(graph), _edge_ref(graph)
+_node_id(graph)
-{
+{}
-// Copy _graph_ref to graph
-_graph.reserveNode(countNodes(_graph_ref) + 1);
-_graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
-copyGraph(_graph, _graph_ref)
-.edgeMap(_capacity, capacity)
-.edgeMap(_cost, cost)
-.nodeRef(_node_ref)
-.edgeRef(_edge_ref)
-.run();
-_supply[_node_ref[s]] =  flow_value;
-_supply[_node_ref[t]] = -flow_value;
-_valid_supply = true;
-}
 /// Destructor.
 ~NetworkSimplex() {
 if (_local_flow) delete _flow_result;
 if (_local_potential) delete _potential_result;
 /// function is \f$ O(e) \f$.
 ///
 /// \pre \ref run() must be called before using this function.
 Cost totalCost() const {
 Cost c = 0;
-for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
+for (EdgeIt e(_orig_graph); e != INVALID; ++e)
-c += (*_flow_result)[e] * _cost[_edge_ref[e]];
+c += (*_flow_result)[e] * _orig_cost[e];
 return c;
 }
 /// @}
 private:
-// Extend the underlying graph and initialize all the node and edge maps
+// Initialize internal data structures
 bool init() {
-if (!_valid_supply) return false;
 // Initialize result maps
 if (!_flow_result) {
-_flow_result = new FlowMap(_graph_ref);
+_flow_result = new FlowMap(_orig_graph);
 _local_flow = true;
 }
 if (!_potential_result) {
-_potential_result = new PotentialMap(_graph_ref);
+_potential_result = new PotentialMap(_orig_graph);
 _local_potential = true;
 }
-// Initialize the edge vector and the edge maps
+// Initialize vectors
-_edges.reserve(countEdges(_graph));
+_node_num = countNodes(_orig_graph);
-for (EdgeIt e(_graph); e != INVALID; ++e) {
+_edge_num = countEdges(_orig_graph);
-_edges.push_back(e);
+int all_node_num = _node_num + 1;
-_flow[e] = 0;
+int all_edge_num = _edge_num + _node_num;
-_state[e] = STATE_LOWER;
-}
+_edge.resize(_edge_num);
+_node.reserve(_node_num);
-// Add an artificial root node to the graph
+_source.resize(all_edge_num);
-_root = _graph.addNode();
+_target.resize(all_edge_num);
-_parent[_root] = INVALID;
-_pred_edge[_root] = INVALID;
+_cap.resize(all_edge_num);
+_cost.resize(all_edge_num);
+_supply.resize(all_node_num);
+_flow.resize(all_edge_num, 0);
+_pi.resize(all_node_num, 0);
+_depth.resize(all_node_num);
+_parent.resize(all_node_num);
+_pred.resize(all_node_num);
+_thread.resize(all_node_num);
+_forward.resize(all_node_num);
+_state.resize(all_edge_num, STATE_LOWER);
+// Initialize node related data
+bool valid_supply = true;
+if (_orig_supply) {
+Supply sum = 0;
+int i = 0;
+for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
+_node.push_back(n);
+_node_id[n] = i;
+_supply[i] = (*_orig_supply)[n];
+sum += _supply[i];
+}
+valid_supply = (sum == 0);
+} else {
+int i = 0;
+for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
+_node.push_back(n);
+_node_id[n] = i;
+_supply[i] = 0;
+}
+_supply[_node_id[_orig_source]] =  _orig_flow_value;
+_supply[_node_id[_orig_target]] = -_orig_flow_value;
+}
+if (!valid_supply) return false;
+// Set data for the artificial root node
+_root = _node_num;
 _depth[_root] = 0;
+_parent[_root] = -1;
+_pred[_root] = -1;
+_thread[_root] = 0;
 _supply[_root] = 0;
-_potential[_root] = 0;
+_pi[_root] = 0;
-// Add artificial edges to the graph and initialize the node maps
+// Store the edges in a mixed order
-// of the spanning tree data structure
+int k = std::max(int(sqrt(_edge_num)), 10);
-Node last = _root;
+int i = 0;
-Edge e;
+for (EdgeIt e(_orig_graph); e != INVALID; ++e) {
+_edge[i] = e;
+if ((i += k) >= _edge_num) i = (i % k) + 1;
+}
+// Initialize edge maps
+for (int i = 0; i != _edge_num; ++i) {
+Edge e = _edge[i];
+_source[i] = _node_id[_orig_graph.source(e)];
+_target[i] = _node_id[_orig_graph.target(e)];
+_cost[i] = _orig_cost[e];
+_cap[i] = _orig_cap[e];
+}
+// Remove non-zero lower bounds
+if (_orig_lower) {
+for (int i = 0; i != _edge_num; ++i) {
+Capacity c = (*_orig_lower)[_edge[i]];
+if (c != 0) {
+_cap[i] -= c;
+_supply[_source[i]] -= c;
+_supply[_target[i]] += c;
+}
+}
+}
+// Add artificial edges and initialize the spanning tree data structure
 Cost max_cost = std::numeric_limits<Cost>::max() / 4;
-for (NodeIt u(_graph); u != INVALID; ++u) {
+Capacity max_cap = std::numeric_limits<Capacity>::max();
-if (u == _root) continue;
+for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) {
-_thread[last] = u;
+_thread[u] = u + 1;
-last = u;
+_depth[u] = 1;
 _parent[u] = _root;
-_depth[u] = 1;
+_pred[u] = e;
 if (_supply[u] >= 0) {
-e = _graph.addEdge(u, _root);
 _flow[e] = _supply[u];
 _forward[u] = true;
-_potential[u] = -max_cost;
+_pi[u] = -max_cost;
 } else {
-e = _graph.addEdge(_root, u);
 _flow[e] = -_supply[u];
 _forward[u] = false;
-_potential[u] = max_cost;
+_pi[u] = max_cost;
 }
 _cost[e] = max_cost;
-_capacity[e] = std::numeric_limits<Capacity>::max();
+_cap[e] = max_cap;
 _state[e] = STATE_TREE;
-_pred_edge[u] = e;
+}
-}
-_thread[last] = _root;
 return true;
 }
 // Find the join node
 void findJoinNode() {
-Node u = _graph.source(_in_edge);
+int u = _source[_in_edge];
-Node v = _graph.target(_in_edge);
+int v = _target[_in_edge];
+while (_depth[u] > _depth[v]) u = _parent[u];
+while (_depth[v] > _depth[u]) v = _parent[v];
 while (u != v) {
-if (_depth[u] == _depth[v]) {
+u = _parent[u];
-u = _parent[u];
+v = _parent[v];
-v = _parent[v];
-}
-else if (_depth[u] > _depth[v]) u = _parent[u];
-else v = _parent[v];
 }
 join = u;
 }
 // Find the leaving edge of the cycle and returns true if the
 // leaving edge is not the same as the entering edge
 bool findLeavingEdge() {
 // Initialize first and second nodes according to the direction
 // of the cycle
 if (_state[_in_edge] == STATE_LOWER) {
-first  = _graph.source(_in_edge);
+first  = _source[_in_edge];
-second = _graph.target(_in_edge);
+second = _target[_in_edge];
 } else {
-first  = _graph.target(_in_edge);
+first  = _target[_in_edge];
-second = _graph.source(_in_edge);
+second = _source[_in_edge];
 }
-delta = _capacity[_in_edge];
+delta = _cap[_in_edge];
-bool result = false;
+int result = 0;
 Capacity d;
-Edge e;
+int e;
 // Search the cycle along the path form the first node to the root
-for (Node u = first; u != join; u = _parent[u]) {
+for (int u = first; u != join; u = _parent[u]) {
-e = _pred_edge[u];
+e = _pred[u];
-d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
+d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
 if (d < delta) {
 delta = d;
 u_out = u;
-u_in = first;
+result = 1;
-v_in = second;
-result = true;
 }
 }
 // Search the cycle along the path form the second node to the root
-for (Node u = second; u != join; u = _parent[u]) {
+for (int u = second; u != join; u = _parent[u]) {
-e = _pred_edge[u];
+e = _pred[u];
-d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
+d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
 if (d <= delta) {
 delta = d;
 u_out = u;
-u_in = second;
+result = 2;
-v_in = first;
+}
-result = true;
+}
-}
-}
+if (result == 1) {
-return result;
+u_in = first;
-}
+v_in = second;
+} else {
-// Change _flow and _state edge maps
+u_in = second;
-void changeFlows(bool change) {
+v_in = first;
+}
+return result != 0;
+}
+// Change _flow and _state vectors
+void changeFlow(bool change) {
 // Augment along the cycle
 if (delta > 0) {
 Capacity val = _state[_in_edge] * delta;
 _flow[_in_edge] += val;
-for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {
+for (int u = _source[_in_edge]; u != join; u = _parent[u]) {
-_flow[_pred_edge[u]] += _forward[u] ? -val : val;
+_flow[_pred[u]] += _forward[u] ? -val : val;
 }
-for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {
+for (int u = _target[_in_edge]; u != join; u = _parent[u]) {
-_flow[_pred_edge[u]] += _forward[u] ? val : -val;
+_flow[_pred[u]] += _forward[u] ? val : -val;
 }
 }
 // Update the state of the entering and leaving edges
 if (change) {
 _state[_in_edge] = STATE_TREE;
-_state[_pred_edge[u_out]] =
+_state[_pred[u_out]] =
-(_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
+(_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
 } else {
 _state[_in_edge] = -_state[_in_edge];
 }
 }
-// Update _thread and _parent node maps
+// Update _thread and _parent vectors
 void updateThreadParent() {
-Node u;
+int u;
 v_out = _parent[u_out];
 // Handle the case when join and v_out coincide
 bool par_first = false;
 if (join == v_out) {
 for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
 _thread[u] = right;
 }
 }
-// Update _pred_edge and _forward node maps.
+// Update _pred and _forward vectors
 void updatePredEdge() {
-Node u = u_out, v;
+int u = u_out, v;
 while (u != u_in) {
 v = _parent[u];
-_pred_edge[u] = _pred_edge[v];
+_pred[u] = _pred[v];
 _forward[u] = !_forward[v];
 u = v;
 }
-_pred_edge[u_in] = _in_edge;
+_pred[u_in] = _in_edge;
-_forward[u_in] = (u_in == _graph.source(_in_edge));
+_forward[u_in] = (u_in == _source[_in_edge]);
 }
-// Update _depth and _potential node maps
+// Update _depth and _potential vectors
 void updateDepthPotential() {
 _depth[u_in] = _depth[v_in] + 1;
 Cost sigma = _forward[u_in] ?
-_potential[v_in] - _potential[u_in] - _cost[_pred_edge[u_in]] :
+_pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
-_potential[v_in] - _potential[u_in] + _cost[_pred_edge[u_in]];
+_pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
-_potential[u_in] += sigma;
+_pi[u_in] += sigma;
-for(Node u = _thread[u_in]; _parent[u] != INVALID; u = _thread[u]) {
+for(int u = _thread[u_in]; _parent[u] != -1; u = _thread[u]) {
 _depth[u] = _depth[_parent[u]] + 1;
 if (_depth[u] <= _depth[u_in]) break;
-_potential[u] += sigma;
+_pi[u] += sigma;
 }
 }
 // Execute the algorithm
 bool start(PivotRuleEnum pivot_rule) {
 return false;
 }
 template<class PivotRuleImplementation>
 bool start() {
-PivotRuleImplementation pivot(*this, _edges);
+PivotRuleImplementation pivot(*this);
 // Execute the network simplex algorithm
 while (pivot.findEnteringEdge()) {
 findJoinNode();
 bool change = findLeavingEdge();
-changeFlows(change);
+changeFlow(change);
 if (change) {
 updateThreadParent();
 updatePredEdge();
 updateDepthPotential();
 }
 }
 // Check if the flow amount equals zero on all the artificial edges
-for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
+for (int e = _edge_num; e != _edge_num + _node_num; ++e) {
 if (_flow[e] > 0) return false;
-for (OutEdgeIt e(_graph, _root); e != INVALID; ++e)
+}
-if (_flow[e] > 0) return false;
 // Copy flow values to _flow_result
-if (_lower) {
+if (_orig_lower) {
-for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
+for (int i = 0; i != _edge_num; ++i) {
-(*_flow_result)[e] = (*_lower)[e] + _flow[_edge_ref[e]];
+Edge e = _edge[i];
+(*_flow_result)[e] = (*_orig_lower)[e] + _flow[i];
+}
 } else {
-for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
+for (int i = 0; i != _edge_num; ++i) {
-(*_flow_result)[e] = _flow[_edge_ref[e]];
+(*_flow_result)[_edge[i]] = _flow[i];
+}
 }
 // Copy potential values to _potential_result
-for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
+for (int i = 0; i != _node_num; ++i) {
-(*_potential_result)[n] = _potential[_node_ref[n]];
+(*_potential_result)[_node[i]] = _pi[i];
+}
 return true;
 }
 }; //class NetworkSimplex

changeset 2634	e98bbe64cca4
parent 2630	d239741cfd44
child 2635	23570e368afa