[Lemon-commits] kpeter: r3458 - lemon/trunk/lemon
Lemon SVN
svn at lemon.cs.elte.hu
Mon Feb 18 04:32:06 CET 2008
Author: kpeter
Date: Mon Feb 18 04:32:06 2008
New Revision: 3458
Modified:
lemon/trunk/lemon/network_simplex.h
Log:
Major improvements in NetworkSimplex.
Main changes:
- Use -potenital[] instead of potential[] to conform to the usual
terminology.
- Use function parameter instead of #define commands to select pivot rule.
- Use much faster implementation for the candidate list pivot rule.
It is about 5-20 times faster now.
- Add a new pivot rule called "Limited Search" that is a modified
version of "Block Search". It is about 25 percent faster on rather
sparse graphs.
- By default "Limited Search" is used for sparse graphs and
"Block Search" is used otherwise. This combined method is the most
efficient on every input class.
- Change the name of private members to start with "_".
- Change the name of function parameters not to start with "_".
- Remove unnecessary documentation for private members.
- Many doc improvements.
Modified: lemon/trunk/lemon/network_simplex.h
==============================================================================
--- lemon/trunk/lemon/network_simplex.h (original)
+++ lemon/trunk/lemon/network_simplex.h Mon Feb 18 04:32:06 2008
@@ -22,47 +22,15 @@
/// \ingroup min_cost_flow
///
/// \file
-/// \brief The network simplex algorithm for finding a minimum cost flow.
+/// \brief Network simplex algorithm for finding a minimum cost flow.
+#include <vector>
#include <limits>
+
#include <lemon/graph_adaptor.h>
#include <lemon/graph_utils.h>
#include <lemon/smart_graph.h>
-
-/// \brief The pivot rule used in the algorithm.
-//#define FIRST_ELIGIBLE_PIVOT
-//#define BEST_ELIGIBLE_PIVOT
-#define EDGE_BLOCK_PIVOT
-//#define CANDIDATE_LIST_PIVOT
-//#define SORTED_LIST_PIVOT
-
-//#define _DEBUG_ITER_
-
-
-// State constant for edges at their lower bounds.
-#define LOWER 1
-// State constant for edges in the spanning tree.
-#define TREE 0
-// State constant for edges at their upper bounds.
-#define UPPER -1
-
-#ifdef EDGE_BLOCK_PIVOT
- #include <lemon/math.h>
- #define MIN_BLOCK_SIZE 10
-#endif
-
-#ifdef CANDIDATE_LIST_PIVOT
- #include <vector>
- #define LIST_LENGTH_DIV 50
- #define MINOR_LIMIT_DIV 200
-#endif
-
-#ifdef SORTED_LIST_PIVOT
- #include <vector>
- #include <algorithm>
- #define LIST_LENGTH_DIV 100
- #define LOWER_DIV 4
-#endif
+#include <lemon/math.h>
namespace lemon {
@@ -75,31 +43,38 @@
/// \ref NetworkSimplex implements the network simplex algorithm for
/// finding a minimum cost flow.
///
- /// \param Graph The directed graph type the algorithm runs on.
- /// \param LowerMap The type of the lower bound map.
- /// \param CapacityMap The type of the capacity (upper bound) map.
- /// \param CostMap The type of the cost (length) map.
- /// \param SupplyMap The type of the supply map.
+ /// \tparam Graph The directed graph type the algorithm runs on.
+ /// \tparam LowerMap The type of the lower bound map.
+ /// \tparam CapacityMap The type of the capacity (upper bound) map.
+ /// \tparam CostMap The type of the cost (length) map.
+ /// \tparam SupplyMap The type of the supply map.
///
/// \warning
- /// - Edge capacities and costs should be non-negative integers.
- /// However \c CostMap::Value should be signed type.
- /// - Supply values should be signed integers.
- /// - \c LowerMap::Value must be convertible to
- /// \c CapacityMap::Value and \c CapacityMap::Value must be
- /// convertible to \c SupplyMap::Value.
+ /// - Edge capacities and costs should be \e non-negative \e integers.
+ /// - Supply values should be \e signed \e integers.
+ /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
+ /// - \c CapacityMap::Value and \c SupplyMap::Value must be
+ /// convertible to each other.
+ /// - All value types must be convertible to \c CostMap::Value, which
+ /// must be signed type.
+ ///
+ /// \note \ref NetworkSimplex provides six different pivot rule
+ /// implementations that significantly affect the efficiency of the
+ /// algorithm.
+ /// By default a combined pivot rule is used, which is the fastest
+ /// implementation according to our benchmark tests.
+ /// Another pivot rule can be selected using \ref run() function
+ /// with the proper parameter.
///
/// \author Peter Kovacs
template < typename Graph,
typename LowerMap = typename Graph::template EdgeMap<int>,
- typename CapacityMap = LowerMap,
+ typename CapacityMap = typename Graph::template EdgeMap<int>,
typename CostMap = typename Graph::template EdgeMap<int>,
- typename SupplyMap = typename Graph::template NodeMap
- <typename CapacityMap::Value> >
+ typename SupplyMap = typename Graph::template NodeMap<int> >
class NetworkSimplex
{
- typedef typename LowerMap::Value Lower;
typedef typename CapacityMap::Value Capacity;
typedef typename CostMap::Value Cost;
typedef typename SupplyMap::Value Supply;
@@ -107,7 +82,6 @@
typedef SmartGraph SGraph;
GRAPH_TYPEDEFS(typename SGraph);
- typedef typename SGraph::template EdgeMap<Lower> SLowerMap;
typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
typedef typename SGraph::template EdgeMap<Cost> SCostMap;
typedef typename SGraph::template NodeMap<Supply> SSupplyMap;
@@ -129,77 +103,380 @@
/// The type of the potential map.
typedef typename Graph::template NodeMap<Cost> PotentialMap;
- protected:
+ public:
+
+ /// Enum type to select the pivot rule used by \ref run().
+ enum PivotRuleEnum {
+ FIRST_ELIGIBLE_PIVOT,
+ BEST_ELIGIBLE_PIVOT,
+ BLOCK_SEARCH_PIVOT,
+ LIMITED_SEARCH_PIVOT,
+ CANDIDATE_LIST_PIVOT,
+ COMBINED_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;
+ const SGraph &_gr;
+ const SCostMap &_cost_map;
+ const SPotentialMap &_pot_map;
public:
- ReducedCostMap( const SGraph &_gr,
- const SCostMap &_cm,
- const SPotentialMap &_pm ) :
- gr(_gr), cost_map(_cm), pot_map(_pm) {}
+ ///\e
+ ReducedCostMap( const SGraph &gr,
+ const SCostMap &cost_map,
+ const SPotentialMap &pot_map ) :
+ _gr(gr), _cost_map(cost_map), _pot_map(pm) {}
+ ///\e
Cost operator[](const Edge &e) const {
- return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
+ return _cost_map[e] + _pot_map[_gr.source(e)]
+ - _pot_map[_gr.target(e)];
}
}; //class ReducedCostMap
- protected:
+ 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:
+
+ NetworkSimplex &_ns;
+ EdgeIt _next_edge;
+
+ public:
+
+ /// Constructor.
+ FirstEligiblePivotRule(NetworkSimplex &ns) :
+ _ns(ns), _next_edge(ns._graph) {}
+
+ /// Finds the next entering edge.
+ bool findEnteringEdge() {
+ for (EdgeIt e = _next_edge; e != INVALID; ++e) {
+ if (_ns._state[e] * _ns._red_cost[e] < 0) {
+ _ns._in_edge = e;
+ _next_edge = ++e;
+ return true;
+ }
+ }
+ for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
+ if (_ns._state[e] * _ns._red_cost[e] < 0) {
+ _ns._in_edge = e;
+ _next_edge = ++e;
+ 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
+ /// for the \ref NetworkSimplex "network simplex" algorithm.
+ class BestEligiblePivotRule
+ {
+ private:
+
+ NetworkSimplex &_ns;
+
+ public:
- /// The directed graph the algorithm runs on.
- SGraph graph;
- /// The original graph.
- const Graph &graph_ref;
- /// The original lower bound map.
- const LowerMap *lower;
- /// The capacity map.
- SCapacityMap capacity;
- /// The cost map.
- SCostMap cost;
- /// The supply map.
- SSupplyMap supply;
- /// The reduced cost map.
- ReducedCostMap red_cost;
- bool valid_supply;
-
- /// The edge map of the current flow.
- SCapacityMap flow;
- /// The potential node map.
- SPotentialMap potential;
- FlowMap flow_result;
- PotentialMap potential_result;
-
- /// Node reference for the original graph.
- NodeRefMap node_ref;
- /// Edge reference for the original graph.
- EdgeRefMap edge_ref;
-
- /// The \c depth node map of the spanning tree structure.
- IntNodeMap depth;
- /// The \c parent node map of the spanning tree structure.
- NodeNodeMap parent;
- /// The \c pred_edge node map of the spanning tree structure.
- EdgeNodeMap pred_edge;
- /// The \c thread node map of the spanning tree structure.
- NodeNodeMap thread;
- /// The \c forward node map of the spanning tree structure.
- BoolNodeMap forward;
- /// The state edge map.
- IntEdgeMap state;
- /// The root node of the starting spanning tree.
- Node root;
+ /// Constructor.
+ BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {}
+
+ /// Finds the next entering edge.
+ bool findEnteringEdge() {
+ Cost min = 0;
+ for (EdgeIt e(_ns._graph); e != INVALID; ++e) {
+ if (_ns._state[e] * _ns._red_cost[e] < min) {
+ min = _ns._state[e] * _ns._red_cost[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
+ /// for the \ref NetworkSimplex "network simplex" algorithm.
+ class BlockSearchPivotRule
+ {
+ private:
+
+ NetworkSimplex &_ns;
+ EdgeIt _next_edge, _min_edge;
+ int _block_size;
+
+ static const int MIN_BLOCK_SIZE = 10;
+
+ public:
+
+ /// Constructor.
+ BlockSearchPivotRule(NetworkSimplex &ns) :
+ _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
+ {
+ _block_size = 2 * int(sqrt(countEdges(_ns._graph)));
+ if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE;
+ }
+
+ /// Finds the next entering edge.
+ bool findEnteringEdge() {
+ Cost curr, min = 0;
+ int cnt = 0;
+ for (EdgeIt e = _next_edge; e != INVALID; ++e) {
+ if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
+ min = curr;
+ _min_edge = e;
+ }
+ if (++cnt == _block_size) {
+ if (min < 0) break;
+ cnt = 0;
+ }
+ }
+ if (min == 0) {
+ for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
+ if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
+ min = curr;
+ _min_edge = e;
+ }
+ if (++cnt == _block_size) {
+ if (min < 0) break;
+ cnt = 0;
+ }
+ }
+ }
+ _ns._in_edge = _min_edge;
+ _next_edge = ++_min_edge;
+ return min < 0;
+ }
+ }; //class BlockSearchPivotRule
+
+ /// \brief Implementation of the "Limited Search" pivot rule for the
+ /// \ref NetworkSimplex "network simplex" algorithm.
+ ///
+ /// This class implements the "Limited Search" pivot rule
+ /// for the \ref NetworkSimplex "network simplex" algorithm.
+ class LimitedSearchPivotRule
+ {
+ private:
+
+ NetworkSimplex &_ns;
+ EdgeIt _next_edge, _min_edge;
+ int _sample_size;
+
+ static const int MIN_SAMPLE_SIZE = 10;
+ static const double SAMPLE_SIZE_FACTOR = 0.0015;
+
+ public:
+
+ /// Constructor.
+ LimitedSearchPivotRule(NetworkSimplex &ns) :
+ _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
+ {
+ _sample_size = int(SAMPLE_SIZE_FACTOR * countEdges(_ns._graph));
+ if (_sample_size < MIN_SAMPLE_SIZE)
+ _sample_size = MIN_SAMPLE_SIZE;
+ }
+
+ /// Finds the next entering edge.
+ bool findEnteringEdge() {
+ Cost curr, min = 0;
+ int cnt = 0;
+ for (EdgeIt e = _next_edge; e != INVALID; ++e) {
+ if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
+ min = curr;
+ _min_edge = e;
+ }
+ if (curr < 0 && ++cnt == _sample_size) break;
+ }
+ if (min == 0) {
+ for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
+ if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
+ min = curr;
+ _min_edge = e;
+ }
+ if (curr < 0 && ++cnt == _sample_size) break;
+ }
+ }
+ _ns._in_edge = _min_edge;
+ _next_edge = ++_min_edge;
+ return min < 0;
+ }
+ }; //class LimitedSearchPivotRule
+
+ /// \brief Implementation of the "Candidate List" pivot rule for the
+ /// \ref NetworkSimplex "network simplex" algorithm.
+ ///
+ /// This class implements the "Candidate List" pivot rule
+ /// for the \ref NetworkSimplex "network simplex" algorithm.
+ class CandidateListPivotRule
+ {
+ private:
+
+ NetworkSimplex &_ns;
+
+ // The list of candidate edges.
+ std::vector<Edge> _candidates;
+ // The maximum length of the edge list.
+ int _list_length;
+ // The maximum number of minor iterations between two major
+ // itarations.
+ int _minor_limit;
+
+ int _minor_count;
+ EdgeIt _next_edge;
+
+ static const double LIST_LENGTH_FACTOR = 0.002;
+ static const double MINOR_LIMIT_FACTOR = 0.1;
+ static const int MIN_LIST_LENGTH = 10;
+ static const int MIN_MINOR_LIMIT = 2;
+
+ public:
+
+ /// Constructor.
+ CandidateListPivotRule(NetworkSimplex &ns) :
+ _ns(ns), _next_edge(ns._graph)
+ {
+ int edge_num = countEdges(_ns._graph);
+ _minor_count = 0;
+ _list_length = int(edge_num * LIST_LENGTH_FACTOR);
+ if (_list_length < MIN_LIST_LENGTH)
+ _list_length = MIN_LIST_LENGTH;
+ _minor_limit = int(_list_length * MINOR_LIMIT_FACTOR);
+ if (_minor_limit < MIN_MINOR_LIMIT)
+ _minor_limit = MIN_MINOR_LIMIT;
+ }
+
+ /// Finds the next entering edge.
+ bool findEnteringEdge() {
+ Cost min, curr;
+ if (_minor_count < _minor_limit && _candidates.size() > 0) {
+ // Minor iteration
+ ++_minor_count;
+ Edge e;
+ min = 0;
+ for (int i = 0; i < int(_candidates.size()); ++i) {
+ e = _candidates[i];
+ if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
+ min = curr;
+ _ns._in_edge = e;
+ }
+ }
+ if (min < 0) return true;
+ }
+
+ // Major iteration
+ _candidates.clear();
+ EdgeIt e = _next_edge;
+ min = 0;
+ for ( ; e != INVALID; ++e) {
+ if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
+ _candidates.push_back(e);
+ if (curr < min) {
+ min = curr;
+ _ns._in_edge = e;
+ }
+ if (int(_candidates.size()) == _list_length) break;
+ }
+ }
+ if (int(_candidates.size()) < _list_length) {
+ for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) {
+ if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
+ _candidates.push_back(e);
+ if (curr < min) {
+ min = curr;
+ _ns._in_edge = e;
+ }
+ if (int(_candidates.size()) == _list_length) break;
+ }
+ }
+ }
+ if (_candidates.size() == 0) return false;
+ _minor_count = 1;
+ _next_edge = ++e;
+ return true;
+ }
+ }; //class CandidateListPivotRule
+
+ private:
+
+ // State constant for edges at their lower bounds.
+ static const int STATE_LOWER = 1;
+ // State constant for edges in the spanning tree.
+ static const int STATE_TREE = 0;
+ // State constant for edges at their upper bounds.
+ static const int STATE_UPPER = -1;
+
+ // Constant for the combined pivot rule.
+ static const int COMBINED_PIVOT_MAX_DEG = 5;
+
+ private:
+
+ // The directed graph the algorithm runs on
+ SGraph _graph;
+ // The original graph
+ const Graph &_graph_ref;
+ // 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 root node of the starting spanning tree
+ Node _root;
+
+ // The reduced cost map
+ ReducedCostMap _red_cost;
+
+ // Members for handling the original graph
+ FlowMap _flow_result;
+ PotentialMap _potential_result;
+ NodeRefMap _node_ref;
+ EdgeRefMap _edge_ref;
// The entering edge of the current pivot iteration.
- Edge in_edge;
+ Edge _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;
@@ -208,82 +485,56 @@
// pivot iteration.
Capacity delta;
-#ifdef EDGE_BLOCK_PIVOT
- /// The size of blocks for the "Edge Block" pivot rule.
- int block_size;
-#endif
-#ifdef CANDIDATE_LIST_PIVOT
- /// \brief The list of candidate edges for the "Candidate List"
- /// pivot rule.
- std::vector<Edge> candidates;
- /// \brief The maximum length of the edge list for the
- /// "Candidate List" pivot rule.
- int list_length;
- /// \brief The maximum number of minor iterations between two major
- /// itarations.
- int minor_limit;
- /// \brief The number of minor iterations.
- int minor_count;
-#endif
-#ifdef SORTED_LIST_PIVOT
- /// \brief The list of candidate edges for the "Sorted List"
- /// pivot rule.
- std::vector<Edge> candidates;
- /// \brief The maximum length of the edge list for the
- /// "Sorted List" pivot rule.
- int list_length;
- int list_index;
-#endif
-
public :
/// \brief General constructor of the class (with lower bounds).
///
/// General constructor of the class (with lower bounds).
///
- /// \param _graph The directed graph the algorithm runs on.
- /// \param _lower The lower bounds of the edges.
- /// \param _capacity The capacities (upper bounds) of the edges.
- /// \param _cost The cost (length) values of the edges.
- /// \param _supply The supply values of the nodes (signed).
- NetworkSimplex( const Graph &_graph,
- const LowerMap &_lower,
- const CapacityMap &_capacity,
- const CostMap &_cost,
- const SupplyMap &_supply ) :
- graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
- supply(graph), flow(graph), flow_result(_graph), potential(graph),
- potential_result(_graph), depth(graph), parent(graph),
- pred_edge(graph), thread(graph), forward(graph), state(graph),
- node_ref(graph_ref), edge_ref(graph_ref),
- red_cost(graph, cost, potential)
+ /// \param graph The directed graph the algorithm runs on.
+ /// \param lower The lower bounds of the edges.
+ /// \param capacity The capacities (upper bounds) of the edges.
+ /// \param cost The cost (length) values of the edges.
+ /// \param supply The supply values of the nodes (signed).
+ NetworkSimplex( const Graph &graph,
+ const LowerMap &lower,
+ 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),
+ _flow_result(graph), _potential_result(graph),
+ _node_ref(graph), _edge_ref(graph)
{
// Checking 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;
-
- // Copying graph_ref to graph
- graph.reserveNode(countNodes(graph_ref) + 1);
- graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref));
- copyGraph(graph, graph_ref)
- .edgeMap(cost, _cost)
- .nodeRef(node_ref)
- .edgeRef(edge_ref)
+ for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
+ sum += supply[n];
+ if (!(_valid_supply = sum == 0)) return;
+
+ // Copying _graph_ref to _graph
+ _graph.reserveNode(countNodes(_graph_ref) + 1);
+ _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
+ copyGraph(_graph, _graph_ref)
+ .edgeMap(_cost, cost)
+ .nodeRef(_node_ref)
+ .edgeRef(_edge_ref)
.run();
// Removing non-zero lower bounds
- for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
- capacity[edge_ref[e]] = _capacity[e] - _lower[e];
+ for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
+ _capacity[_edge_ref[e]] = capacity[e] - lower[e];
}
- for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
- Supply s = _supply[n];
- for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
- s += _lower[e];
- for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
- s -= _lower[e];
- supply[node_ref[n]] = s;
+ for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
+ Supply s = supply[n];
+ for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
+ s += lower[e];
+ for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
+ s -= lower[e];
+ _supply[_node_ref[n]] = s;
}
}
@@ -291,34 +542,35 @@
///
/// General constructor of the class (without lower bounds).
///
- /// \param _graph The directed graph the algorithm runs on.
- /// \param _capacity The capacities (upper bounds) of the edges.
- /// \param _cost The cost (length) values of the edges.
- /// \param _supply The supply values of the nodes (signed).
- NetworkSimplex( const Graph &_graph,
- const CapacityMap &_capacity,
- const CostMap &_cost,
- const SupplyMap &_supply ) :
- graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
- supply(graph), flow(graph), flow_result(_graph), potential(graph),
- potential_result(_graph), depth(graph), parent(graph),
- pred_edge(graph), thread(graph), forward(graph), state(graph),
- node_ref(graph_ref), edge_ref(graph_ref),
- red_cost(graph, cost, potential)
+ /// \param graph The directed graph the algorithm runs on.
+ /// \param capacity The capacities (upper bounds) of the edges.
+ /// \param cost The cost (length) values of the edges.
+ /// \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),
+ _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),
+ _flow_result(graph), _potential_result(graph),
+ _node_ref(graph), _edge_ref(graph)
{
// Checking 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;
-
- // Copying graph_ref to graph
- copyGraph(graph, graph_ref)
- .edgeMap(capacity, _capacity)
- .edgeMap(cost, _cost)
- .nodeMap(supply, _supply)
- .nodeRef(node_ref)
- .edgeRef(edge_ref)
+ for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
+ sum += supply[n];
+ if (!(_valid_supply = sum == 0)) return;
+
+ // Copying _graph_ref to graph
+ copyGraph(_graph, _graph_ref)
+ .edgeMap(_capacity, capacity)
+ .edgeMap(_cost, cost)
+ .nodeMap(_supply, supply)
+ .nodeRef(_node_ref)
+ .edgeRef(_edge_ref)
.run();
}
@@ -326,115 +578,150 @@
///
/// Simple constructor of the class (with lower bounds).
///
- /// \param _graph The directed graph the algorithm runs on.
- /// \param _lower The lower bounds of the edges.
- /// \param _capacity The capacities (upper bounds) of the edges.
- /// \param _cost The cost (length) values of the edges.
- /// \param _s The source node.
- /// \param _t The target node.
- /// \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 LowerMap &_lower,
- const CapacityMap &_capacity,
- const CostMap &_cost,
- typename Graph::Node _s,
- typename Graph::Node _t,
- typename SupplyMap::Value _flow_value ) :
- graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
- supply(graph), flow(graph), flow_result(_graph), potential(graph),
- potential_result(_graph), depth(graph), parent(graph),
- pred_edge(graph), thread(graph), forward(graph), state(graph),
- node_ref(graph_ref), edge_ref(graph_ref),
- red_cost(graph, cost, potential)
+ /// \param graph The directed graph the algorithm runs on.
+ /// \param lower The lower bounds of the edges.
+ /// \param capacity The capacities (upper bounds) of the edges.
+ /// \param cost The cost (length) values of the edges.
+ /// \param s The source node.
+ /// \param t The target node.
+ /// \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 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), _flow(_graph),
+ _potential(_graph), _depth(_graph), _parent(_graph),
+ _pred_edge(_graph), _thread(_graph), _forward(_graph),
+ _state(_graph), _red_cost(_graph, _cost, _potential),
+ _flow_result(graph), _potential_result(graph),
+ _node_ref(graph), _edge_ref(graph)
{
- // Copying graph_ref to graph
- copyGraph(graph, graph_ref)
- .edgeMap(cost, _cost)
- .nodeRef(node_ref)
- .edgeRef(edge_ref)
+ // Copying _graph_ref to graph
+ copyGraph(_graph, _graph_ref)
+ .edgeMap(_cost, cost)
+ .nodeRef(_node_ref)
+ .edgeRef(_edge_ref)
.run();
// Removing non-zero lower bounds
- for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
- capacity[edge_ref[e]] = _capacity[e] - _lower[e];
+ for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
+ _capacity[_edge_ref[e]] = capacity[e] - lower[e];
}
- for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
- Supply s = 0;
- if (n == _s) s = _flow_value;
- if (n == _t) s = -_flow_value;
- for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
- s += _lower[e];
- for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
- s -= _lower[e];
- supply[node_ref[n]] = s;
+ for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
+ Supply sum = 0;
+ if (n == s) sum = flow_value;
+ if (n == t) sum = -flow_value;
+ for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
+ sum += lower[e];
+ for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
+ sum -= lower[e];
+ _supply[_node_ref[n]] = sum;
}
- valid_supply = true;
+ _valid_supply = true;
}
/// \brief Simple constructor of the class (without lower bounds).
///
/// Simple constructor of the class (without lower bounds).
///
- /// \param _graph The directed graph the algorithm runs on.
- /// \param _capacity The capacities (upper bounds) of the edges.
- /// \param _cost The cost (length) values of the edges.
- /// \param _s The source node.
- /// \param _t The target node.
- /// \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,
- typename Graph::Node _t,
- typename SupplyMap::Value _flow_value ) :
- graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
- supply(graph, 0), flow(graph), flow_result(_graph), potential(graph),
- potential_result(_graph), depth(graph), parent(graph),
- pred_edge(graph), thread(graph), forward(graph), state(graph),
- node_ref(graph_ref), edge_ref(graph_ref),
- red_cost(graph, cost, potential)
+ /// \param graph The directed graph the algorithm runs on.
+ /// \param capacity The capacities (upper bounds) of the edges.
+ /// \param cost The cost (length) values of the edges.
+ /// \param s The source node.
+ /// \param t The target node.
+ /// \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,
+ 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),
+ _flow_result(graph), _potential_result(graph),
+ _node_ref(graph), _edge_ref(graph)
{
- // Copying graph_ref to graph
- copyGraph(graph, graph_ref)
- .edgeMap(capacity, _capacity)
- .edgeMap(cost, _cost)
- .nodeRef(node_ref)
- .edgeRef(edge_ref)
+ // Copying _graph_ref to graph
+ 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;
+ _supply[_node_ref[s]] = flow_value;
+ _supply[_node_ref[t]] = -flow_value;
+ _valid_supply = true;
}
/// \brief Runs the algorithm.
///
/// Runs the algorithm.
///
+ /// \param pivot_rule The pivot rule that is used during the
+ /// algorithm.
+ ///
+ /// The available pivot rules:
+ ///
+ /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in
+ /// a wraparound fashion in every iteration
+ /// (\ref FirstEligiblePivotRule).
+ ///
+ /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in
+ /// every iteration (\ref BestEligiblePivotRule).
+ ///
+ /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in
+ /// every iteration in a wraparound fashion and the best eligible
+ /// edge is selected from this block (\ref BlockSearchPivotRule).
+ ///
+ /// - LIMITED_SEARCH_PIVOT A specified number of eligible edges are
+ /// examined in every iteration in a wraparound fashion and the best
+ /// one is selected from them (\ref LimitedSearchPivotRule).
+ ///
+ /// - CANDIDATE_LIST_PIVOT In major iterations a candidate list is
+ /// built from eligible edges and it is used for edge selection in
+ /// the following minor iterations (\ref CandidateListPivotRule).
+ ///
+ /// - COMBINED_PIVOT This is a combined version of the two fastest
+ /// pivot rules.
+ /// For rather sparse graphs \ref LimitedSearchPivotRule
+ /// "Limited Search" implementation is used, otherwise
+ /// \ref BlockSearchPivotRule "Block Search" pivot rule is used.
+ /// According to our benchmark tests this combined method is the
+ /// most efficient.
+ ///
/// \return \c true if a feasible flow can be found.
- bool run() {
- return init() && start();
+ bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) {
+ return init() && start(pivot_rule);
}
- /// \brief Returns a const reference to the flow map.
+ /// \brief Returns a const reference to the edge map storing the
+ /// found flow.
///
- /// Returns a const reference to the flow map.
+ /// Returns a const reference to the edge map storing the found flow.
///
/// \pre \ref run() must be called before using this function.
const FlowMap& flowMap() const {
- return flow_result;
+ return _flow_result;
}
- /// \brief Returns a const reference to the potential map (the dual
- /// solution).
+ /// \brief Returns a const reference to the node map storing the
+ /// found potentials (the dual solution).
///
- /// Returns a const reference to the potential map (the dual
- /// solution).
+ /// Returns a const reference to the node map storing the found
+ /// potentials (the dual solution).
///
/// \pre \ref run() must be called before using this function.
const PotentialMap& potentialMap() const {
- return potential_result;
+ return _potential_result;
}
/// \brief Returns the total cost of the found flow.
@@ -445,248 +732,75 @@
/// \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 (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
+ c += _flow_result[e] * _cost[_edge_ref[e]];
return c;
}
- protected:
+ private:
/// \brief Extends the underlaying graph and initializes all the
/// node and edge maps.
bool init() {
- if (!valid_supply) return false;
+ if (!_valid_supply) return false;
// Initializing state and flow maps
- for (EdgeIt e(graph); e != INVALID; ++e) {
- flow[e] = 0;
- state[e] = LOWER;
+ for (EdgeIt e(_graph); e != INVALID; ++e) {
+ _flow[e] = 0;
+ _state[e] = STATE_LOWER;
}
// Adding an artificial root node to the graph
- root = graph.addNode();
- parent[root] = INVALID;
- pred_edge[root] = INVALID;
- depth[root] = 0;
- supply[root] = 0;
- potential[root] = 0;
+ _root = _graph.addNode();
+ _parent[_root] = INVALID;
+ _pred_edge[_root] = INVALID;
+ _depth[_root] = 0;
+ _supply[_root] = 0;
+ _potential[_root] = 0;
// Adding artificial edges to the graph and initializing the node
// maps of the spanning tree data structure
- Supply sum = 0;
- Node last = root;
+ Node last = _root;
Edge e;
Cost max_cost = std::numeric_limits<Cost>::max() / 4;
- for (NodeIt u(graph); u != INVALID; ++u) {
- if (u == root) continue;
- thread[last] = u;
+ for (NodeIt u(_graph); u != INVALID; ++u) {
+ if (u == _root) continue;
+ _thread[last] = u;
last = u;
- parent[u] = root;
- depth[u] = 1;
- sum += supply[u];
- if (supply[u] >= 0) {
- e = graph.addEdge(u, root);
- flow[e] = supply[u];
- forward[u] = true;
- potential[u] = max_cost;
+ _parent[u] = _root;
+ _depth[u] = 1;
+ if (_supply[u] >= 0) {
+ e = _graph.addEdge(u, _root);
+ _flow[e] = _supply[u];
+ _forward[u] = true;
+ _potential[u] = -max_cost;
} else {
- e = graph.addEdge(root, u);
- flow[e] = -supply[u];
- forward[u] = false;
- potential[u] = -max_cost;
- }
- cost[e] = max_cost;
- capacity[e] = std::numeric_limits<Capacity>::max();
- state[e] = TREE;
- pred_edge[u] = e;
- }
- thread[last] = root;
-
-#ifdef EDGE_BLOCK_PIVOT
- // Initializing block_size for the edge block pivot rule
- int edge_num = countEdges(graph);
- block_size = 2 * int(sqrt(countEdges(graph)));
- if (block_size < MIN_BLOCK_SIZE) block_size = MIN_BLOCK_SIZE;
-#endif
-#ifdef CANDIDATE_LIST_PIVOT
- int edge_num = countEdges(graph);
- minor_count = 0;
- list_length = edge_num / LIST_LENGTH_DIV;
- minor_limit = edge_num / MINOR_LIMIT_DIV;
-#endif
-#ifdef SORTED_LIST_PIVOT
- int edge_num = countEdges(graph);
- list_index = 0;
- list_length = edge_num / LIST_LENGTH_DIV;
-#endif
-
- return sum == 0;
- }
-
-#ifdef FIRST_ELIGIBLE_PIVOT
- /// \brief Finds entering edge according to the "First Eligible"
- /// pivot rule.
- bool findEnteringEdge(EdgeIt &next_edge) {
- for (EdgeIt e = next_edge; e != INVALID; ++e) {
- if (state[e] * red_cost[e] < 0) {
- in_edge = e;
- next_edge = ++e;
- return true;
- }
- }
- for (EdgeIt e(graph); e != next_edge; ++e) {
- if (state[e] * red_cost[e] < 0) {
- in_edge = e;
- next_edge = ++e;
- return true;
- }
- }
- return false;
- }
-#endif
-
-#ifdef BEST_ELIGIBLE_PIVOT
- /// \brief Finds entering edge according to the "Best Eligible"
- /// pivot rule.
- bool findEnteringEdge() {
- Cost min = 0;
- for (EdgeIt e(graph); e != INVALID; ++e) {
- if (state[e] * red_cost[e] < min) {
- min = state[e] * red_cost[e];
- in_edge = e;
- }
- }
- return min < 0;
- }
-#endif
-
-#ifdef EDGE_BLOCK_PIVOT
- /// \brief Finds entering edge according to the "Edge Block"
- /// pivot rule.
- bool findEnteringEdge(EdgeIt &next_edge) {
- // Performing edge block selection
- Cost curr, min = 0;
- EdgeIt min_edge(graph);
- int cnt = 0;
- for (EdgeIt e = next_edge; e != INVALID; ++e) {
- if ((curr = state[e] * red_cost[e]) < min) {
- min = curr;
- min_edge = e;
- }
- if (++cnt == block_size) {
- if (min < 0) break;
- cnt = 0;
- }
- }
- if (!(min < 0)) {
- for (EdgeIt e(graph); e != next_edge; ++e) {
- if ((curr = state[e] * red_cost[e]) < min) {
- min = curr;
- min_edge = e;
- }
- if (++cnt == block_size) {
- if (min < 0) break;
- cnt = 0;
- }
- }
- }
- in_edge = min_edge;
- if ((next_edge = ++min_edge) == INVALID)
- next_edge = EdgeIt(graph);
- return min < 0;
- }
-#endif
-
-#ifdef CANDIDATE_LIST_PIVOT
- /// \brief Finds entering edge according to the "Candidate List"
- /// pivot rule.
- bool findEnteringEdge() {
- typedef typename std::vector<Edge>::iterator ListIt;
-
- if (minor_count >= minor_limit || candidates.size() == 0) {
- // Major iteration
- candidates.clear();
- for (EdgeIt e(graph); e != INVALID; ++e) {
- if (state[e] * red_cost[e] < 0) {
- candidates.push_back(e);
- if (candidates.size() == list_length) break;
- }
- }
- if (candidates.size() == 0) return false;
+ e = _graph.addEdge(_root, u);
+ _flow[e] = -_supply[u];
+ _forward[u] = false;
+ _potential[u] = max_cost;
+ }
+ _cost[e] = max_cost;
+ _capacity[e] = std::numeric_limits<Capacity>::max();
+ _state[e] = STATE_TREE;
+ _pred_edge[u] = e;
}
+ _thread[last] = _root;
- // Minor iteration
- ++minor_count;
- Cost min = 0;
- Edge e;
- for (int i = 0; i < candidates.size(); ++i) {
- e = candidates[i];
- if (state[e] * red_cost[e] < min) {
- min = state[e] * red_cost[e];
- in_edge = e;
- }
- }
return true;
}
-#endif
-#ifdef SORTED_LIST_PIVOT
- /// \brief Functor class to compare edges during sort of the
- /// candidate list.
- class SortFunc
- {
- private:
- const IntEdgeMap &st;
- const ReducedCostMap &rc;
- public:
- SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) :
- st(_st), rc(_rc) {}
- bool operator()(const Edge &e1, const Edge &e2) {
- return st[e1] * rc[e1] < st[e2] * rc[e2];
- }
- };
-
- /// \brief Finds entering edge according to the "Sorted List"
- /// pivot rule.
- bool findEnteringEdge() {
- static SortFunc sort_func(state, red_cost);
-
- // Minor iteration
- while (list_index < candidates.size()) {
- in_edge = candidates[list_index++];
- if (state[in_edge] * red_cost[in_edge] < 0) return true;
- }
-
- // Major iteration
- candidates.clear();
- Cost curr, min = 0;
- for (EdgeIt e(graph); e != INVALID; ++e) {
- if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) {
- candidates.push_back(e);
- if (curr < min) min = curr;
- if (candidates.size() == list_length) break;
- }
- }
- if (candidates.size() == 0) return false;
- sort(candidates.begin(), candidates.end(), sort_func);
- in_edge = candidates[0];
- list_index = 1;
- return true;
- }
-#endif
-
- /// \brief Finds the join node.
+ /// Finds the join node.
Node findJoinNode() {
- // Finding the join node
- Node u = graph.source(in_edge);
- Node v = graph.target(in_edge);
+ Node u = _graph.source(_in_edge);
+ Node v = _graph.target(_in_edge);
while (u != v) {
- if (depth[u] == depth[v]) {
- u = parent[u];
- v = parent[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];
+ else if (_depth[u] > _depth[v]) u = _parent[u];
+ else v = _parent[v];
}
return u;
}
@@ -696,23 +810,23 @@
bool findLeavingEdge() {
// Initializing first and second nodes according to the direction
// of the cycle
- if (state[in_edge] == LOWER) {
- first = graph.source(in_edge);
- second = graph.target(in_edge);
+ if (_state[_in_edge] == STATE_LOWER) {
+ first = _graph.source(_in_edge);
+ second = _graph.target(_in_edge);
} else {
- first = graph.target(in_edge);
- second = graph.source(in_edge);
+ first = _graph.target(_in_edge);
+ second = _graph.source(_in_edge);
}
- delta = capacity[in_edge];
+ delta = _capacity[_in_edge];
bool result = false;
Capacity d;
Edge e;
// Searching the cycle along the path form the first node to the
// root node
- for (Node u = first; u != join; u = parent[u]) {
- e = pred_edge[u];
- d = forward[u] ? flow[e] : capacity[e] - flow[e];
+ for (Node u = first; u != join; u = _parent[u]) {
+ e = _pred_edge[u];
+ d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
if (d < delta) {
delta = d;
u_out = u;
@@ -723,9 +837,9 @@
}
// Searching the cycle along the path form the second node to the
// root node
- for (Node u = second; u != join; u = parent[u]) {
- e = pred_edge[u];
- d = forward[u] ? capacity[e] - flow[e] : flow[e];
+ for (Node u = second; u != join; u = _parent[u]) {
+ e = _pred_edge[u];
+ d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
if (d <= delta) {
delta = d;
u_out = u;
@@ -737,134 +851,150 @@
return result;
}
- /// \brief Changes \c flow and \c state edge maps.
+ /// Changes \c flow and \c state edge maps.
void changeFlows(bool change) {
// Augmenting 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;
+ 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 (Node u = graph.target(in_edge); u != join; u = parent[u]) {
- flow[pred_edge[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;
}
}
// Updating the state of the entering and leaving edges
if (change) {
- state[in_edge] = TREE;
- state[pred_edge[u_out]] =
- (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER;
+ _state[_in_edge] = STATE_TREE;
+ _state[_pred_edge[u_out]] =
+ (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
} else {
- state[in_edge] = -state[in_edge];
+ _state[_in_edge] = -_state[_in_edge];
}
}
- /// \brief Updates \c thread and \c parent node maps.
+ /// Updates \c thread and \c parent node maps.
void updateThreadParent() {
Node u;
- v_out = parent[u_out];
+ v_out = _parent[u_out];
// Handling the case when join and v_out coincide
bool par_first = false;
if (join == v_out) {
- for (u = join; u != u_in && u != v_in; u = thread[u]) ;
+ for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
if (u == v_in) {
par_first = true;
- while (thread[u] != u_out) u = thread[u];
+ while (_thread[u] != u_out) u = _thread[u];
first = u;
}
}
// Finding the last successor of u_in (u) and the node after it
// (right) according to the thread index
- for (u = u_in; depth[thread[u]] > depth[u_in]; u = thread[u]) ;
- right = thread[u];
- if (thread[v_in] == u_out) {
- for (last = u; depth[last] > depth[u_out]; last = thread[last]) ;
- if (last == u_out) last = thread[last];
+ for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
+ right = _thread[u];
+ if (_thread[v_in] == u_out) {
+ for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
+ if (last == u_out) last = _thread[last];
}
- else last = thread[v_in];
+ else last = _thread[v_in];
// Updating stem nodes
- thread[v_in] = stem = u_in;
+ _thread[v_in] = stem = u_in;
par_stem = v_in;
while (stem != u_out) {
- thread[u] = new_stem = parent[stem];
+ _thread[u] = new_stem = _parent[stem];
// Finding the node just before the stem node (u) according to
// the original thread index
- for (u = new_stem; thread[u] != stem; u = thread[u]) ;
- thread[u] = right;
+ for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
+ _thread[u] = right;
// Changing the parent node of stem and shifting stem and
// par_stem nodes
- parent[stem] = par_stem;
+ _parent[stem] = par_stem;
par_stem = stem;
stem = new_stem;
// Finding the last successor of stem (u) and the node after it
// (right) according to the thread index
- for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ;
- right = thread[u];
+ for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
+ right = _thread[u];
}
- parent[u_out] = par_stem;
- thread[u] = last;
+ _parent[u_out] = par_stem;
+ _thread[u] = last;
if (join == v_out && par_first) {
- if (first != v_in) thread[first] = right;
+ if (first != v_in) _thread[first] = right;
} else {
- for (u = v_out; thread[u] != u_out; u = thread[u]) ;
- thread[u] = right;
+ for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
+ _thread[u] = right;
}
}
- /// \brief Updates \c pred_edge and \c forward node maps.
+ /// Updates \c pred_edge and \c forward node maps.
void updatePredEdge() {
Node u = u_out, v;
while (u != u_in) {
- v = parent[u];
- pred_edge[u] = pred_edge[v];
- forward[u] = !forward[v];
+ v = _parent[u];
+ _pred_edge[u] = _pred_edge[v];
+ _forward[u] = !_forward[v];
u = v;
}
- pred_edge[u_in] = in_edge;
- forward[u_in] = (u_in == graph.source(in_edge));
+ _pred_edge[u_in] = _in_edge;
+ _forward[u_in] = (u_in == _graph.source(_in_edge));
}
- /// \brief Updates \c depth and \c potential node maps.
+ /// Updates \c depth and \c potential node maps.
void updateDepthPotential() {
- depth[u_in] = depth[v_in] + 1;
- potential[u_in] = forward[u_in] ?
- potential[v_in] + cost[pred_edge[u_in]] :
- potential[v_in] - cost[pred_edge[u_in]];
+ _depth[u_in] = _depth[v_in] + 1;
+ _potential[u_in] = _forward[u_in] ?
+ _potential[v_in] - _cost[_pred_edge[u_in]] :
+ _potential[v_in] + _cost[_pred_edge[u_in]];
- Node u = thread[u_in], v;
+ Node u = _thread[u_in], v;
while (true) {
- v = parent[u];
+ v = _parent[u];
if (v == INVALID) break;
- depth[u] = depth[v] + 1;
- potential[u] = forward[u] ?
- potential[v] + cost[pred_edge[u]] :
- potential[v] - cost[pred_edge[u]];
- if (depth[u] <= depth[v_in]) break;
- u = thread[u];
+ _depth[u] = _depth[v] + 1;
+ _potential[u] = _forward[u] ?
+ _potential[v] - _cost[_pred_edge[u]] :
+ _potential[v] + _cost[_pred_edge[u]];
+ if (_depth[u] <= _depth[v_in]) break;
+ u = _thread[u];
}
}
- /// \brief Executes the algorithm.
+ /// Executes the algorithm.
+ bool start(PivotRuleEnum pivot_rule) {
+ switch (pivot_rule) {
+ case FIRST_ELIGIBLE_PIVOT:
+ return start<FirstEligiblePivotRule>();
+ case BEST_ELIGIBLE_PIVOT:
+ return start<BestEligiblePivotRule>();
+ case BLOCK_SEARCH_PIVOT:
+ return start<BlockSearchPivotRule>();
+ case LIMITED_SEARCH_PIVOT:
+ return start<LimitedSearchPivotRule>();
+ case CANDIDATE_LIST_PIVOT:
+ return start<CandidateListPivotRule>();
+ case COMBINED_PIVOT:
+ if ( countEdges(_graph) / countNodes(_graph) <=
+ COMBINED_PIVOT_MAX_DEG )
+ return start<LimitedSearchPivotRule>();
+ else
+ return start<BlockSearchPivotRule>();
+ }
+ return false;
+ }
+
+ template<class PivotRuleImplementation>
bool start() {
- // Processing pivots
-#ifdef _DEBUG_ITER_
- int iter_num = 0;
-#endif
-#if defined(FIRST_ELIGIBLE_PIVOT) || defined(EDGE_BLOCK_PIVOT)
- EdgeIt next_edge(graph);
- while (findEnteringEdge(next_edge))
-#else
- while (findEnteringEdge())
-#endif
- {
+ PivotRuleImplementation pivot(*this);
+
+ // Executing the network simplex algorithm
+ while (pivot.findEnteringEdge()) {
join = findJoinNode();
bool change = findLeavingEdge();
changeFlows(change);
@@ -873,34 +1003,26 @@
updatePredEdge();
updateDepthPotential();
}
-#ifdef _DEBUG_ITER_
- ++iter_num;
-#endif
- }
-
-#ifdef _DEBUG_ITER_
- std::cout << "Network Simplex algorithm finished. " << iter_num
- << " pivot iterations performed." << std::endl;
-#endif
-
- // Checking 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)
- if (flow[e] > 0) return false;
-
- // Copying 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]];
+ }
+
+ // Checking 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)
+ if (_flow[e] > 0) return false;
+
+ // Copying 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]];
} else {
- for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
- flow_result[e] = flow[edge_ref[e]];
+ for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
+ _flow_result[e] = _flow[_edge_ref[e]];
}
- // Copying potential values to potential_result
- for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
- potential_result[n] = potential[node_ref[n]];
+ // Copying potential values to _potential_result
+ for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
+ _potential_result[n] = _potential[_node_ref[n]];
return true;
}
More information about the Lemon-commits
mailing list