[Lemon-commits] kpeter: r3520 - lemon/trunk/lemon
Lemon SVN
svn at lemon.cs.elte.hu
Fri Feb 6 22:52:34 CET 2009
Author: kpeter
Date: Fri Feb 6 22:52:34 2009
New Revision: 3520
Modified:
lemon/trunk/lemon/network_simplex.h
Log:
Rework Network Simplex
Use simpler and faster graph implementation instead of SmartGraph
Modified: lemon/trunk/lemon/network_simplex.h
==============================================================================
--- lemon/trunk/lemon/network_simplex.h (original)
+++ lemon/trunk/lemon/network_simplex.h Fri Feb 6 22:52:34 2009
@@ -28,9 +28,7 @@
#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 {
@@ -72,29 +70,21 @@
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:
@@ -116,35 +106,6 @@
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.
///
@@ -157,40 +118,52 @@
private:
// References to the NetworkSimplex class
- NetworkSimplex &_ns;
- EdgeVector &_edges;
+ const EdgeVector &_edge;
+ 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) :
- _ns(ns), _edges(edges), _next_edge(0) {}
+ FirstEligiblePivotRule(NetworkSimplex &ns) :
+ _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;
- for (int i = _next_edge; i < int(_edges.size()); ++i) {
- e = _edges[i];
- if (_ns._state[e] * _ns._red_cost[e] < 0) {
- _ns._in_edge = e;
- _next_edge = i + 1;
+ Cost c;
+ for (int e = _next_edge; e < _edge_num; ++e) {
+ c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+ if (c < 0) {
+ _in_edge = e;
+ _next_edge = e + 1;
return true;
}
}
- for (int i = 0; i < _next_edge; ++i) {
- e = _edges[i];
- if (_ns._state[e] * _ns._red_cost[e] < 0) {
- _ns._in_edge = e;
- _next_edge = i + 1;
+ for (int e = 0; e < _next_edge; ++e) {
+ c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+ if (c < 0) {
+ _in_edge = e;
+ _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.
///
@@ -203,30 +176,40 @@
private:
// References to the NetworkSimplex class
- NetworkSimplex &_ns;
- EdgeVector &_edges;
+ const EdgeVector &_edge;
+ 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) :
- _ns(ns), _edges(edges) {}
+ BestEligiblePivotRule(NetworkSimplex &ns) :
+ _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;
- Edge e;
- for (int i = 0; i < int(_edges.size()); ++i) {
- e = _edges[i];
- if (_ns._state[e] * _ns._red_cost[e] < min) {
- min = _ns._state[e] * _ns._red_cost[e];
- _ns._in_edge = e;
+ Cost c, min = 0;
+ for (int e = 0; e < _edge_num; ++e) {
+ c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+ if (c < min) {
+ min = c;
+ _in_edge = e;
}
}
return min < 0;
}
+
}; //class BestEligiblePivotRule
+
/// \brief Implementation of the "Block Search" pivot rule for the
/// \ref NetworkSimplex "network simplex" algorithm.
///
@@ -239,37 +222,45 @@
private:
// References to the NetworkSimplex class
- NetworkSimplex &_ns;
- EdgeVector &_edges;
+ const EdgeVector &_edge;
+ 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) :
- _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
+ BlockSearchPivotRule(NetworkSimplex &ns) :
+ _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;
- Edge e;
+ Cost c, min = 0;
int cnt = _block_size;
- int i;
- for (i = _next_edge; i < int(_edges.size()); ++i) {
- e = _edges[i];
- if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
- min = curr;
- _min_edge = i;
+ int e, min_edge = _next_edge;
+ for (e = _next_edge; e < _edge_num; ++e) {
+ c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+ if (c < min) {
+ min = c;
+ min_edge = e;
}
if (--cnt == 0) {
if (min < 0) break;
@@ -277,11 +268,11 @@
}
}
if (min == 0 || cnt > 0) {
- for (i = 0; i < _next_edge; ++i) {
- e = _edges[i];
- if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
- min = curr;
- _min_edge = i;
+ for (e = 0; e < _next_edge; ++e) {
+ c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+ if (c < min) {
+ min = c;
+ min_edge = e;
}
if (--cnt == 0) {
if (min < 0) break;
@@ -290,12 +281,14 @@
}
}
if (min >= 0) return false;
- _ns._in_edge = _edges[_min_edge];
- _next_edge = i;
+ _in_edge = min_edge;
+ _next_edge = e;
return true;
}
+
}; //class BlockSearchPivotRule
+
/// \brief Implementation of the "Candidate List" pivot rule for the
/// \ref NetworkSimplex "network simplex" algorithm.
///
@@ -308,19 +301,28 @@
private:
// References to the NetworkSimplex class
- NetworkSimplex &_ns;
- EdgeVector &_edges;
+ const EdgeVector &_edge;
+ const IntVector &_source;
+ const IntVector &_target;
+ const CostVector &_cost;
+ const IntVector &_state;
+ const CostVector &_pi;
+ int &_in_edge;
+ int _edge_num;
- EdgeVector _candidates;
+ // 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) :
- _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
+ CandidateListPivotRule(NetworkSimplex &ns) :
+ _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;
@@ -328,7 +330,7 @@
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 );
@@ -338,51 +340,52 @@
/// 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];
- if (curr < min) {
- min = curr;
- _ns._in_edge = e;
+ c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+ if (c < min) {
+ min = c;
+ 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 (i = _next_edge; i < int(_edges.size()); ++i) {
- e = _edges[i];
- if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
+ for (e = _next_edge; e < _edge_num; ++e) {
+ c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+ if (c < 0) {
_candidates[_curr_length++] = e;
- if (curr < min) {
- min = curr;
- _min_edge = i;
+ if (c < min) {
+ min = c;
+ min_edge = e;
}
if (_curr_length == _list_length) break;
}
}
if (_curr_length < _list_length) {
- for (i = 0; i < _next_edge; ++i) {
- e = _edges[i];
- if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
+ for (e = 0; e < _next_edge; ++e) {
+ c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+ if (c < 0) {
_candidates[_curr_length++] = e;
- if (curr < min) {
- min = curr;
- _min_edge = i;
+ if (c < min) {
+ min = c;
+ min_edge = e;
}
if (_curr_length == _list_length) break;
}
@@ -390,12 +393,14 @@
}
if (_curr_length == 0) return false;
_minor_count = 1;
- _ns._in_edge = _edges[_min_edge];
- _next_edge = i;
+ _in_edge = min_edge;
+ _next_edge = e;
return true;
}
+
}; //class CandidateListPivotRule
+
/// \brief Implementation of the "Altering Candidate List" pivot rule
/// for the \ref NetworkSimplex "network simplex" algorithm.
///
@@ -408,23 +413,29 @@
private:
// References to the NetworkSimplex class
- NetworkSimplex &_ns;
- EdgeVector &_edges;
+ const EdgeVector &_edge;
+ const IntVector &_source;
+ const IntVector &_target;
+ const CostVector &_cost;
+ const IntVector &_state;
+ const CostVector &_pi;
+ int &_in_edge;
+ int _edge_num;
- EdgeVector _candidates;
- SCostMap _cand_cost;
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) {}
- bool operator()(const Edge &e1, const Edge &e2) {
- return _map[e1] > _map[e2];
+ SortFunc(const CostVector &map) : _map(map) {}
+ bool operator()(int left, int right) {
+ return _map[left] > _map[right];
}
};
@@ -433,9 +444,11 @@
public:
/// Constructor
- AlteringListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
- _ns(ns), _edges(edges), _cand_cost(_ns._graph),
- _next_edge(0), _sort_func(_cand_cost)
+ AlteringListPivotRule(NetworkSimplex &ns) :
+ _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), _cand_cost(ns._edge_num),_sort_func(_cand_cost)
{
// The main parameters of the pivot rule
const double BLOCK_SIZE_FACTOR = 1.5;
@@ -443,7 +456,7 @@
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 );
@@ -454,11 +467,13 @@
/// Find next entering edge
bool findEnteringEdge() {
// Check the current candidate list
- Edge e;
- for (int ix = 0; ix < _curr_length; ++ix) {
- e = _candidates[ix];
- if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) {
- _candidates[ix--] = _candidates[--_curr_length];
+ int e;
+ for (int i = 0; i < _curr_length; ++i) {
+ e = _candidates[i];
+ _cand_cost[e] = _state[e] *
+ (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+ if (_cand_cost[e] >= 0) {
+ _candidates[i--] = _candidates[--_curr_length];
}
}
@@ -466,11 +481,13 @@
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];
- if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
+
+ for (int e = _next_edge; e < _edge_num; ++e) {
+ _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;
@@ -479,11 +496,12 @@
}
}
if (_curr_length <= limit) {
- for (int i = 0; i < _next_edge; ++i) {
- e = _edges[i];
- if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
+ for (int e = 0; e < _next_edge; ++e) {
+ _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;
@@ -500,12 +518,13 @@
_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:
@@ -519,66 +538,87 @@
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;
- // The capacity map
- SCapacityMap _capacity;
- // The cost map
- SCostMap _cost;
- // The supply map
- SSupplyMap _supply;
- bool _valid_supply;
-
- // Edge map of the current flow
- SCapacityMap _flow;
- // Node map of the current potentials
- SPotentialMap _potential;
-
- // The depth node map of the spanning tree structure
- 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;
+ const LowerMap *_orig_lower;
+ // The original capacity map
+ const CapacityMap &_orig_cap;
+ // The original cost map
+ const CostMap &_orig_cost;
+ // The original supply map
+ const SupplyMap *_orig_supply;
+ // The source node (if no supply map was given)
+ Node _orig_source;
+ // The target node (if no supply map was given)
+ Node _orig_target;
+ // The flow value (if no supply map was given)
+ Capacity _orig_flow_value;
- // The non-artifical edges
- EdgeVector _edges;
-
- // Members for handling the original graph
+ // The flow result map
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;
- Node right, first, second, last;
- Node stem, par_stem, new_stem;
+ int join, u_in, v_in, u_out, v_out;
+ int right, first, second, last;
+ 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).
///
@@ -594,41 +634,12 @@
const CapacityMap &capacity,
const CostMap &cost,
const SupplyMap &supply ) :
- _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
- _cost(_graph), _supply(_graph), _flow(_graph),
- _potential(_graph), _depth(_graph), _parent(_graph),
- _pred_edge(_graph), _thread(_graph), _forward(_graph),
- _state(_graph), _red_cost(_graph, _cost, _potential),
+ _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
+ _orig_cost(cost), _orig_supply(&supply),
_flow_result(NULL), _potential_result(NULL),
_local_flow(false), _local_potential(false),
- _node_ref(graph), _edge_ref(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];
- }
- }
- }
+ _node_id(graph)
+ {}
/// \brief General constructor (without lower bounds).
///
@@ -642,32 +653,12 @@
const CapacityMap &capacity,
const CostMap &cost,
const SupplyMap &supply ) :
- _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
- _cost(_graph), _supply(_graph), _flow(_graph),
- _potential(_graph), _depth(_graph), _parent(_graph),
- _pred_edge(_graph), _thread(_graph), _forward(_graph),
- _state(_graph), _red_cost(_graph, _cost, _potential),
+ _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
+ _orig_cost(cost), _orig_supply(&supply),
_flow_result(NULL), _potential_result(NULL),
_local_flow(false), _local_potential(false),
- _node_ref(graph), _edge_ref(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();
- }
+ _node_id(graph)
+ {}
/// \brief Simple constructor (with lower bounds).
///
@@ -685,40 +676,15 @@
const LowerMap &lower,
const CapacityMap &capacity,
const CostMap &cost,
- typename Graph::Node s,
- typename Graph::Node t,
- typename SupplyMap::Value flow_value ) :
- _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
- _cost(_graph), _supply(_graph, 0), _flow(_graph),
- _potential(_graph), _depth(_graph), _parent(_graph),
- _pred_edge(_graph), _thread(_graph), _forward(_graph),
- _state(_graph), _red_cost(_graph, _cost, _potential),
+ Node s, Node t,
+ Capacity flow_value ) :
+ _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
+ _orig_cost(cost), _orig_supply(NULL),
+ _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
_flow_result(NULL), _potential_result(NULL),
_local_flow(false), _local_potential(false),
- _node_ref(graph), _edge_ref(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;
- }
+ _node_id(graph)
+ {}
/// \brief Simple constructor (without lower bounds).
///
@@ -734,31 +700,15 @@
NetworkSimplex( const Graph &graph,
const CapacityMap &capacity,
const CostMap &cost,
- typename Graph::Node s,
- typename Graph::Node t,
- typename SupplyMap::Value flow_value ) :
- _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
- _cost(_graph), _supply(_graph, 0), _flow(_graph),
- _potential(_graph), _depth(_graph), _parent(_graph),
- _pred_edge(_graph), _thread(_graph), _forward(_graph),
- _state(_graph), _red_cost(_graph, _cost, _potential),
+ Node s, Node t,
+ Capacity flow_value ) :
+ _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
+ _orig_cost(cost), _orig_supply(NULL),
+ _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
_flow_result(NULL), _potential_result(NULL),
_local_flow(false), _local_potential(false),
- _node_ref(graph), _edge_ref(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;
- }
+ _node_id(graph)
+ {}
/// Destructor.
~NetworkSimplex() {
@@ -894,8 +844,8 @@
/// \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)
- c += (*_flow_result)[e] * _cost[_edge_ref[e]];
+ for (EdgeIt e(_orig_graph); e != INVALID; ++e)
+ c += (*_flow_result)[e] * _orig_cost[e];
return c;
}
@@ -903,153 +853,216 @@
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
- _edges.reserve(countEdges(_graph));
- for (EdgeIt e(_graph); e != INVALID; ++e) {
- _edges.push_back(e);
- _flow[e] = 0;
- _state[e] = STATE_LOWER;
+
+ // Initialize vectors
+ _node_num = countNodes(_orig_graph);
+ _edge_num = countEdges(_orig_graph);
+ int all_node_num = _node_num + 1;
+ int all_edge_num = _edge_num + _node_num;
+
+ _edge.resize(_edge_num);
+ _node.reserve(_node_num);
+ _source.resize(all_edge_num);
+ _target.resize(all_edge_num);
+
+ _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;
}
-
- // Add an artificial root node to the graph
- _root = _graph.addNode();
- _parent[_root] = INVALID;
- _pred_edge[_root] = INVALID;
+ 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;
+
+ // Store the edges in a mixed order
+ int k = std::max(int(sqrt(_edge_num)), 10);
+ int i = 0;
+ 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];
+ }
- // Add artificial edges to the graph and initialize the node maps
- // of the spanning tree data structure
- Node last = _root;
- Edge 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) {
- if (u == _root) continue;
- _thread[last] = u;
- last = u;
- _parent[u] = _root;
+ Capacity max_cap = std::numeric_limits<Capacity>::max();
+ for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) {
+ _thread[u] = u + 1;
_depth[u] = 1;
+ _parent[u] = _root;
+ _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);
- Node v = _graph.target(_in_edge);
+ int u = _source[_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];
- v = _parent[v];
- }
- else if (_depth[u] > _depth[v]) u = _parent[u];
- else v = _parent[v];
+ u = _parent[u];
+ v = _parent[v];
}
join = u;
}
- // Find the leaving edge of the cycle and returns true if the
+ // 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);
- second = _graph.target(_in_edge);
+ first = _source[_in_edge];
+ second = _target[_in_edge];
} else {
- first = _graph.target(_in_edge);
- second = _graph.source(_in_edge);
+ first = _target[_in_edge];
+ second = _source[_in_edge];
}
- delta = _capacity[_in_edge];
- bool result = false;
+ delta = _cap[_in_edge];
+ 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]) {
- e = _pred_edge[u];
- d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
+ for (int u = first; u != join; u = _parent[u]) {
+ e = _pred[u];
+ d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
if (d < delta) {
delta = d;
u_out = u;
- u_in = first;
- v_in = second;
- result = true;
+ result = 1;
}
}
// Search the cycle along the path form the second node to the root
- for (Node u = second; u != join; u = _parent[u]) {
- e = _pred_edge[u];
- d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
+ for (int u = second; u != join; u = _parent[u]) {
+ e = _pred[u];
+ d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
if (d <= delta) {
delta = d;
u_out = u;
- u_in = second;
- v_in = first;
- result = true;
+ result = 2;
}
}
- return result;
+
+ if (result == 1) {
+ u_in = first;
+ v_in = second;
+ } else {
+ u_in = second;
+ v_in = first;
+ }
+ return result != 0;
}
- // Change _flow and _state edge maps
- void changeFlows(bool change) {
+ // 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]) {
- _flow[_pred_edge[u]] += _forward[u] ? -val : val;
+ for (int u = _source[_in_edge]; u != join; u = _parent[u]) {
+ _flow[_pred[u]] += _forward[u] ? -val : val;
}
- for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {
- _flow[_pred_edge[u]] += _forward[u] ? val : -val;
+ for (int u = _target[_in_edge]; u != join; u = _parent[u]) {
+ _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]] =
- (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
+ _state[_pred[u_out]] =
+ (_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
@@ -1105,30 +1118,30 @@
}
}
- // 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;
- _forward[u_in] = (u_in == _graph.source(_in_edge));
+ _pred[u_in] = _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]] :
- _potential[v_in] - _potential[u_in] + _cost[_pred_edge[u_in]];
- _potential[u_in] += sigma;
- for(Node u = _thread[u_in]; _parent[u] != INVALID; u = _thread[u]) {
+ _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
+ _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
+ _pi[u_in] += sigma;
+ 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;
}
}
@@ -1152,13 +1165,13 @@
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();
@@ -1167,22 +1180,25 @@
}
// Check if the flow amount equals zero on all the artificial edges
- for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
- if (_flow[e] > 0) return false;
- for (OutEdgeIt e(_graph, _root); e != INVALID; ++e)
+ for (int e = _edge_num; e != _edge_num + _node_num; ++e) {
if (_flow[e] > 0) return false;
+ }
// Copy flow values to _flow_result
- if (_lower) {
- for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
- (*_flow_result)[e] = (*_lower)[e] + _flow[_edge_ref[e]];
+ if (_orig_lower) {
+ for (int i = 0; i != _edge_num; ++i) {
+ Edge e = _edge[i];
+ (*_flow_result)[e] = (*_orig_lower)[e] + _flow[i];
+ }
} else {
- for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
- (*_flow_result)[e] = _flow[_edge_ref[e]];
+ for (int i = 0; i != _edge_num; ++i) {
+ (*_flow_result)[_edge[i]] = _flow[i];
+ }
}
// Copy potential values to _potential_result
- for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
- (*_potential_result)[n] = _potential[_node_ref[n]];
+ for (int i = 0; i != _node_num; ++i) {
+ (*_potential_result)[_node[i]] = _pi[i];
+ }
return true;
}
More information about the Lemon-commits
mailing list