# Changeset 2619:30fb4d68b0e8 in lemon-0.x

Ignore:
Timestamp:
10/05/08 15:36:43 (11 years ago)
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3504
Message:

Improve network simplex algorithm

• Remove "Limited Search" and "Combined" pivot rules.
• Add a new pivot rule "Altering Candidate List".
• Make the edge selection faster in every pivot rule.
• Set the default rule to "Block Search".
• Doc improvements.

The algorithm became about 15-35 percent faster on various input files.
"Block Search" pivot rule proved to be by far the fastest on all inputs.

File:
1 edited

### Legend:

Unmodified
 r2593 #include #define _DEBUG_ namespace lemon { /// @{ /// \brief Implementation of the network simplex algorithm for /// finding a minimum cost flow. /// \brief Implementation of the primal network simplex algorithm /// for finding a minimum cost flow. /// /// \ref NetworkSimplex implements the network simplex algorithm for /// finding a minimum cost flow. /// \ref NetworkSimplex implements the primal network simplex algorithm /// for finding a minimum cost flow. /// /// \tparam Graph The directed graph type the algorithm runs on. /// - \c CostMap::Value must be signed type. /// /// \note \ref NetworkSimplex provides six different pivot rule /// \note \ref NetworkSimplex provides five 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. /// By default "Block Search" pivot rule is used, which proved to be /// by far the most efficient according to our benchmark tests. /// However 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, typedef typename SGraph::template NodeMap EdgeNodeMap; typedef typename SGraph::template EdgeMap IntEdgeMap; typedef typename SGraph::template EdgeMap BoolEdgeMap; typedef typename Graph::template NodeMap NodeRefMap; typedef typename Graph::template EdgeMap EdgeRefMap; typedef std::vector EdgeVector; public: BEST_ELIGIBLE_PIVOT, BLOCK_SEARCH_PIVOT, LIMITED_SEARCH_PIVOT, CANDIDATE_LIST_PIVOT, COMBINED_PIVOT ALTERING_LIST_PIVOT }; /// This class implements the "First Eligible" pivot rule /// for the \ref NetworkSimplex "network simplex" algorithm. /// /// For more information see \ref NetworkSimplex::run(). class FirstEligiblePivotRule { private: // References to the NetworkSimplex class NetworkSimplex &_ns; EdgeIt _next_edge; EdgeVector &_edges; int _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) { /// Constructor FirstEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) : _ns(ns), _edges(edges), _next_edge(0) {} /// Find next entering edge inline 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 = ++e; _next_edge = i + 1; return true; } } for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { 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 = ++e; _next_edge = i + 1; return true; } /// This class implements the "Best Eligible" pivot rule /// for the \ref NetworkSimplex "network simplex" algorithm. /// /// For more information see \ref NetworkSimplex::run(). class BestEligiblePivotRule { private: // References to the NetworkSimplex class NetworkSimplex &_ns; EdgeVector &_edges; public: /// Constructor. BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {} /// Finds the next entering edge. bool findEnteringEdge() { /// Constructor BestEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) : _ns(ns), _edges(edges) {} /// Find next entering edge inline bool findEnteringEdge() { Cost min = 0; for (EdgeIt e(_ns._graph); e != INVALID; ++e) { 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]; /// This class implements the "Block Search" pivot rule /// for the \ref NetworkSimplex "network simplex" algorithm. /// /// For more information see \ref NetworkSimplex::run(). class BlockSearchPivotRule { private: // References to the NetworkSimplex class NetworkSimplex &_ns; EdgeIt _next_edge, _min_edge; EdgeVector &_edges; int _block_size; static const int MIN_BLOCK_SIZE = 10; int _next_edge, _min_edge; public: /// Constructor. BlockSearchPivotRule(NetworkSimplex &ns) : _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph) /// Constructor BlockSearchPivotRule(NetworkSimplex &ns, EdgeVector &edges) : _ns(ns), _edges(edges), _next_edge(0), _min_edge(0) { _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() { // 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())), MIN_BLOCK_SIZE ); } /// Find next entering edge inline bool findEnteringEdge() { Cost curr, min = 0; int cnt = 0; for (EdgeIt e = _next_edge; e != INVALID; ++e) { Edge e; 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 = e; _min_edge = i; } if (++cnt == _block_size) { if (--cnt == 0) { if (min < 0) break; cnt = 0; cnt = _block_size; } } if (min == 0) { for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { 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 = e; _min_edge = i; } if (++cnt == _block_size) { if (--cnt == 0) { if (min < 0) break; cnt = 0; cnt = _block_size; } } } _ns._in_edge = _min_edge; _next_edge = ++_min_edge; return min < 0; if (min >= 0) return false; _ns._in_edge = _edges[_min_edge]; _next_edge = i; return true; } }; //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 SAMPLE_SIZE_FACTOR = 15; static const int MIN_SAMPLE_SIZE = 10; public: /// Constructor. LimitedSearchPivotRule(NetworkSimplex &ns) : _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph) { _sample_size = countEdges(_ns._graph) * SAMPLE_SIZE_FACTOR / 10000; 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 /// This class implements the "Candidate List" pivot rule /// for the \ref NetworkSimplex "network simplex" algorithm. /// /// For more information see \ref NetworkSimplex::run(). class CandidateListPivotRule { private: // References to the NetworkSimplex class NetworkSimplex &_ns; // The list of candidate edges. std::vector _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 int LIST_LENGTH_FACTOR = 20; static const int MINOR_LIMIT_FACTOR = 10; static const int MIN_LIST_LENGTH = 10; static const int MIN_MINOR_LIMIT = 2; EdgeVector &_edges; EdgeVector _candidates; int _list_length, _minor_limit; int _curr_length, _minor_count; int _next_edge, _min_edge; public: /// Constructor. CandidateListPivotRule(NetworkSimplex &ns) : _ns(ns), _next_edge(ns._graph) /// Constructor CandidateListPivotRule(NetworkSimplex &ns, EdgeVector &edges) : _ns(ns), _edges(edges), _next_edge(0), _min_edge(0) { int edge_num = countEdges(_ns._graph); _minor_count = 0; _list_length = edge_num * LIST_LENGTH_FACTOR / 10000; if (_list_length < MIN_LIST_LENGTH) _list_length = MIN_LIST_LENGTH; _minor_limit = _list_length * MINOR_LIMIT_FACTOR / 100; if (_minor_limit < MIN_MINOR_LIMIT) _minor_limit = MIN_MINOR_LIMIT; } /// Finds the next entering edge. bool findEnteringEdge() { // 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())), 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 inline bool findEnteringEdge() { Cost min, curr; if (_minor_count < _minor_limit && _candidates.size() > 0) { // Minor iteration if (_curr_length > 0 && _minor_count < _minor_limit) { // Minor iteration: selecting the best eligible edge from // the current candidate list ++_minor_count; Edge e; min = 0; for (int i = 0; i < int(_candidates.size()); ++i) { for (int i = 0; i < _curr_length; ++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); curr = _ns._state[e] * _ns._red_cost[e]; if (curr < min) { min = curr; _ns._in_edge = e; } if (int(_candidates.size()) == _list_length) break; if (curr >= 0) { _candidates[i--] = _candidates[--_curr_length]; } } } if (int(_candidates.size()) < _list_length) { for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) { if (min < 0) return true; } // Major iteration: building 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) { _candidates[_curr_length++] = e; if (curr < min) { min = curr; _min_edge = i; } 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) { _candidates.push_back(e); _candidates[_curr_length++] = e; if (curr < min) { min = curr; _ns._in_edge = e; _min_edge = i; } if (int(_candidates.size()) == _list_length) break; if (_curr_length == _list_length) break; } } } if (_candidates.size() == 0) return false; if (_curr_length == 0) return false; _minor_count = 1; _next_edge = ++e; _ns._in_edge = _edges[_min_edge]; _next_edge = i; 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 /// for the \ref NetworkSimplex "network simplex" algorithm. /// /// For more information see \ref NetworkSimplex::run(). class AlteringListPivotRule { private: // References to the NetworkSimplex class NetworkSimplex &_ns; EdgeVector &_edges; EdgeVector _candidates; SCostMap _cand_cost; int _block_size, _head_length, _curr_length; int _next_edge; // Functor class to compare edges during sort of the candidate list class SortFunc { private: const SCostMap &_map; public: SortFunc(const SCostMap &map) : _map(map) {} bool operator()(const Edge &e1, const Edge &e2) { return _map[e1] < _map[e2]; } }; SortFunc _sort_func; public: /// Constructor AlteringListPivotRule(NetworkSimplex &ns, EdgeVector &edges) : _ns(ns), _edges(edges), _cand_cost(_ns._graph), _next_edge(0), _sort_func(_cand_cost) { // The main parameters of the pivot rule const double BLOCK_SIZE_FACTOR = 1.0; const int MIN_BLOCK_SIZE = 10; const double HEAD_LENGTH_FACTOR = 0.1; const int MIN_HEAD_LENGTH = 5; _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())), 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 inline bool findEnteringEdge() { // Checking the current candidate list Edge e; for (int idx = 0; idx < _curr_length; ++idx) { e = _candidates[idx]; if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) { _candidates[idx--] = _candidates[--_curr_length]; } } // Extending 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]; if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) { _candidates[_curr_length++] = e; last_edge = i; } 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) { e = _edges[i]; if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) { _candidates[_curr_length++] = e; last_edge = i; } if (--cnt == 0) { if (_curr_length > limit) break; limit = 0; cnt = _block_size; } } } if (_curr_length == 0) return false; _next_edge = last_edge + 1; // Sorting the list partially EdgeVector::iterator sort_end = _candidates.begin(); EdgeVector::iterator vector_end = _candidates.begin(); for (int idx = 0; idx < _curr_length; ++idx) { ++vector_end; if (idx <= _head_length) ++sort_end; } partial_sort(_candidates.begin(), sort_end, vector_end, _sort_func); _ns._in_edge = _candidates[0]; if (_curr_length > _head_length) { _candidates[0] = _candidates[_head_length - 1]; _curr_length = _head_length - 1; } else { _candidates[0] = _candidates[_curr_length - 1]; --_curr_length; } return true; } }; //class AlteringListPivotRule private: STATE_LOWER =  1 }; // Constant for the combined pivot rule. static const int COMBINED_PIVOT_MAX_DEG = 5; private: // The reduced cost map ReducedCostMap _red_cost; // The non-artifical edges EdgeVector _edges; // Members for handling the original graph } /// \brief Sets the flow map. /// /// Sets the flow map. /// \brief Set the flow map. /// /// Set the flow map. /// /// \return \c (*this) } /// \brief Sets the potential map. /// /// Sets the potential map. /// \brief Set the potential map. /// /// Set the potential map. /// /// \return \c (*this) /// \name Execution control /// The only way to execute the algorithm is to call the run() /// function. /// @{ /// 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. /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is /// built from eligible edges in a wraparound fashion and in the /// following minor iterations the best eligible edge is selected /// from this list (\ref CandidateListPivotRule). /// /// - ALTERING_LIST_PIVOT It is a modified version of the /// "Candidate List" pivot rule. It keeps only the several best /// eligible edges from the former candidate list and extends this /// list in every iteration (\ref AlteringListPivotRule). /// /// According to our comprehensive benchmark tests the "Block Search" /// pivot rule proved to be by far the fastest and the most robust /// on various test inputs. Thus it is the default option. /// /// \return \c true if a feasible flow can be found. bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) { bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) { return init() && start(pivot_rule); } /// \name Query Functions /// The result of the algorithm can be obtained using these /// functions. /// \n run() must be called before using them. /// The results of the algorithm can be obtained using these /// functions.\n /// \ref lemon::NetworkSimplex::run() "run()" must be called before /// using them. /// @{ /// \brief Returns a const reference to the edge map storing the /// \brief Return a const reference to the edge map storing the /// found flow. /// /// Returns a const reference to the edge map storing the found flow. /// Return a const reference to the edge map storing the found flow. /// /// \pre \ref run() must be called before using this function. } /// \brief Returns a const reference to the node map storing the /// \brief Return a const reference to the node map storing the /// found potentials (the dual solution). /// /// Returns a const reference to the node map storing the found /// Return a const reference to the node map storing the found /// potentials (the dual solution). /// } /// \brief Returns the flow on the given edge. /// /// Returns the flow on the given edge. /// \brief Return the flow on the given edge. /// /// Return the flow on the given edge. /// /// \pre \ref run() must be called before using this function. } /// \brief Returns the potential of the given node. /// /// Returns the potential of the given node. /// \brief Return the potential of the given node. /// /// Return the potential of the given node. /// /// \pre \ref run() must be called before using this function. } /// \brief Returns the total cost of the found flow. /// /// Returns the total cost of the found flow. The complexity of the /// \brief Return the total cost of the found flow. /// /// Return the total cost of the found flow. The complexity of the /// function is \f$O(e) \f$. /// private: /// \brief Extends the underlaying graph and initializes all the /// \brief Extend the underlying graph and initialize all the /// node and edge maps. bool init() { } // Initializing state and flow maps // Initializing 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; } /// Finds the join node. Node findJoinNode() { /// Find the join node. inline Node findJoinNode() { Node u = _graph.source(_in_edge); Node v = _graph.target(_in_edge); } /// \brief Finds the leaving edge of the cycle. Returns \c true if /// the leaving edge is not the same as the entering edge. bool findLeavingEdge() { /// \brief Find the leaving edge of the cycle. /// \return \c true if the leaving edge is not the same as the /// entering edge. inline bool findLeavingEdge() { // Initializing first and second nodes according to the direction // of the cycle } /// Changes \c flow and \c state edge maps. void changeFlows(bool change) { /// Change \c flow and \c state edge maps. inline void changeFlows(bool change) { // Augmenting along the cycle if (delta > 0) { } /// Updates \c thread and \c parent node maps. void updateThreadParent() { /// Update \c thread and \c parent node maps. inline void updateThreadParent() { Node u; v_out = _parent[u_out]; } /// Updates \c pred_edge and \c forward node maps. void updatePredEdge() { /// Update \c pred_edge and \c forward node maps. inline void updatePredEdge() { Node u = u_out, v; while (u != u_in) { } /// Updates \c depth and \c potential node maps. void updateDepthPotential() { /// Update \c depth and \c potential node maps. inline void updateDepthPotential() { _depth[u_in] = _depth[v_in] + 1; _potential[u_in] = _forward[u_in] ? } /// Executes the algorithm. /// Execute the algorithm. bool start(PivotRuleEnum pivot_rule) { // Selecting the pivot rule implementation switch (pivot_rule) { case FIRST_ELIGIBLE_PIVOT: case BLOCK_SEARCH_PIVOT: return start(); case LIMITED_SEARCH_PIVOT: return start(); case CANDIDATE_LIST_PIVOT: return start(); case COMBINED_PIVOT: if ( countEdges(_graph) / countNodes(_graph) <= COMBINED_PIVOT_MAX_DEG ) return start(); else return start(); case ALTERING_LIST_PIVOT: return start(); } return false; template bool start() { PivotRuleImplementation pivot(*this); PivotRuleImplementation pivot(*this, _edges); // Executing the network simplex algorithm