diff --git a/lemon/network_simplex.h b/lemon/network_simplex.h --- a/lemon/network_simplex.h +++ b/lemon/network_simplex.h @@ -40,7 +40,9 @@ /// for finding a \ref min_cost_flow "minimum cost flow". /// /// \ref NetworkSimplex implements the primal Network Simplex algorithm - /// for finding a \ref min_cost_flow "minimum cost flow". + /// for finding a \ref min_cost_flow "minimum cost flow" + /// \ref amo93networkflows, \ref dantzig63linearprog, + /// \ref kellyoneill91netsimplex. /// This algorithm is a specialized version of the linear programming /// simplex method directly for the minimum cost flow problem. /// It is one of the most efficient solution methods. @@ -48,7 +50,7 @@ /// In general this class is the fastest implementation available /// in LEMON for the minimum cost flow problem. /// Moreover it supports both directions of the supply/demand inequality - /// constraints. For more information see \ref SupplyType. + /// constraints. For more information, see \ref SupplyType. /// /// Most of the parameters of the problem (except for the digraph) /// can be given using separate functions, and the algorithm can be @@ -57,16 +59,16 @@ /// /// \tparam GR The digraph type the algorithm runs on. /// \tparam V The value type used for flow amounts, capacity bounds - /// and supply values in the algorithm. By default it is \c int. + /// and supply values in the algorithm. By default, it is \c int. /// \tparam C The value type used for costs and potentials in the - /// algorithm. By default it is the same as \c V. + /// algorithm. By default, it is the same as \c V. /// /// \warning Both value types must be signed and all input data must /// be integer. /// /// \note %NetworkSimplex provides five different pivot rule /// implementations, from which the most efficient one is used - /// by default. For more information see \ref PivotRule. + /// by default. For more information, see \ref PivotRule. template class NetworkSimplex { @@ -122,35 +124,35 @@ /// \ref NetworkSimplex provides five different pivot rule /// implementations that significantly affect the running time /// of the algorithm. - /// By default \ref BLOCK_SEARCH "Block Search" is used, which + /// By default, \ref BLOCK_SEARCH "Block Search" is used, which /// proved to be the most efficient and the most robust on various /// test inputs according to our benchmark tests. - /// However another pivot rule can be selected using the \ref run() + /// However, another pivot rule can be selected using the \ref run() /// function with the proper parameter. enum PivotRule { - /// The First Eligible pivot rule. + /// The \e First \e Eligible pivot rule. /// The next eligible arc is selected in a wraparound fashion /// in every iteration. FIRST_ELIGIBLE, - /// The Best Eligible pivot rule. + /// The \e Best \e Eligible pivot rule. /// The best eligible arc is selected in every iteration. BEST_ELIGIBLE, - /// The Block Search pivot rule. + /// The \e Block \e Search pivot rule. /// A specified number of arcs are examined in every iteration /// in a wraparound fashion and the best eligible arc is selected /// from this block. BLOCK_SEARCH, - /// The Candidate List pivot rule. + /// The \e Candidate \e List pivot rule. /// In a major iteration a candidate list is built from eligible arcs /// in a wraparound fashion and in the following minor iterations /// the best eligible arc is selected from this list. CANDIDATE_LIST, - /// The Altering Candidate List pivot rule. + /// The \e Altering \e Candidate \e List pivot rule. /// It is a modified version of the Candidate List method. /// It keeps only the several best eligible arcs from the former /// candidate list and extends this list in every iteration. @@ -161,8 +163,6 @@ TEMPLATE_DIGRAPH_TYPEDEFS(GR); - typedef std::vector ArcVector; - typedef std::vector NodeVector; typedef std::vector IntVector; typedef std::vector BoolVector; typedef std::vector ValueVector; @@ -364,33 +364,32 @@ bool findEnteringArc() { Cost c, min = 0; int cnt = _block_size; - int e, min_arc = _next_arc; + int e; for (e = _next_arc; e < _search_arc_num; ++e) { c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < min) { min = c; - min_arc = e; + _in_arc = e; } if (--cnt == 0) { - if (min < 0) break; + if (min < 0) goto search_end; cnt = _block_size; } } - if (min == 0 || cnt > 0) { - for (e = 0; e < _next_arc; ++e) { - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); - if (c < min) { - min = c; - min_arc = e; - } - if (--cnt == 0) { - if (min < 0) break; - cnt = _block_size; - } + for (e = 0; e < _next_arc; ++e) { + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); + if (c < min) { + min = c; + _in_arc = e; + } + if (--cnt == 0) { + if (min < 0) goto search_end; + cnt = _block_size; } } if (min >= 0) return false; - _in_arc = min_arc; + + search_end: _next_arc = e; return true; } @@ -428,7 +427,7 @@ _next_arc(0) { // The main parameters of the pivot rule - const double LIST_LENGTH_FACTOR = 1.0; + const double LIST_LENGTH_FACTOR = 0.25; const int MIN_LIST_LENGTH = 10; const double MINOR_LIMIT_FACTOR = 0.1; const int MIN_MINOR_LIMIT = 3; @@ -445,7 +444,7 @@ /// Find next entering arc bool findEnteringArc() { Cost min, c; - int e, min_arc = _next_arc; + int e; if (_curr_length > 0 && _minor_count < _minor_limit) { // Minor iteration: select the best eligible arc from the // current candidate list @@ -456,16 +455,13 @@ c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < min) { min = c; - min_arc = e; + _in_arc = e; } - if (c >= 0) { + else if (c >= 0) { _candidates[i--] = _candidates[--_curr_length]; } } - if (min < 0) { - _in_arc = min_arc; - return true; - } + if (min < 0) return true; } // Major iteration: build a new candidate list @@ -477,27 +473,26 @@ _candidates[_curr_length++] = e; if (c < min) { min = c; - min_arc = e; + _in_arc = e; } - if (_curr_length == _list_length) break; + if (_curr_length == _list_length) goto search_end; } } - if (_curr_length < _list_length) { - for (e = 0; e < _next_arc; ++e) { - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); - if (c < 0) { - _candidates[_curr_length++] = e; - if (c < min) { - min = c; - min_arc = e; - } - if (_curr_length == _list_length) break; + for (e = 0; e < _next_arc; ++e) { + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); + if (c < 0) { + _candidates[_curr_length++] = e; + if (c < min) { + min = c; + _in_arc = e; } + if (_curr_length == _list_length) goto search_end; } } if (_curr_length == 0) return false; + + search_end: _minor_count = 1; - _in_arc = min_arc; _next_arc = e; return true; } @@ -549,7 +544,7 @@ _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost) { // The main parameters of the pivot rule - const double BLOCK_SIZE_FACTOR = 1.5; + 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 = 3; @@ -578,39 +573,35 @@ // Extend the list int cnt = _block_size; - int last_arc = 0; int limit = _head_length; - for (int e = _next_arc; e < _search_arc_num; ++e) { + for (e = _next_arc; e < _search_arc_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_arc = e; } if (--cnt == 0) { - if (_curr_length > limit) break; + if (_curr_length > limit) goto search_end; limit = 0; cnt = _block_size; } } - if (_curr_length <= limit) { - for (int e = 0; e < _next_arc; ++e) { - _cand_cost[e] = _state[e] * - (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); - if (_cand_cost[e] < 0) { - _candidates[_curr_length++] = e; - last_arc = e; - } - if (--cnt == 0) { - if (_curr_length > limit) break; - limit = 0; - cnt = _block_size; - } + for (e = 0; e < _next_arc; ++e) { + _cand_cost[e] = _state[e] * + (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); + if (_cand_cost[e] < 0) { + _candidates[_curr_length++] = e; + } + if (--cnt == 0) { + if (_curr_length > limit) goto search_end; + limit = 0; + cnt = _block_size; } } if (_curr_length == 0) return false; - _next_arc = last_arc + 1; + + search_end: // Make heap of the candidate list (approximating a partial sort) make_heap( _candidates.begin(), _candidates.begin() + _curr_length, @@ -618,6 +609,7 @@ // Pop the first element of the heap _in_arc = _candidates[0]; + _next_arc = e; pop_heap( _candidates.begin(), _candidates.begin() + _curr_length, _sort_func ); _curr_length = std::min(_head_length, _curr_length - 1); @@ -633,7 +625,11 @@ /// The constructor of the class. /// /// \param graph The digraph the algorithm runs on. - NetworkSimplex(const GR& graph) : + /// \param arc_mixing Indicate if the arcs have to be stored in a + /// mixed order in the internal data structure. + /// In special cases, it could lead to better overall performance, + /// but it is usually slower. Therefore it is disabled by default. + NetworkSimplex(const GR& graph, bool arc_mixing = false) : _graph(graph), _node_id(graph), _arc_id(graph), INF(std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : @@ -671,31 +667,33 @@ _last_succ.resize(all_node_num); _state.resize(max_arc_num); - // Copy the graph (store the arcs in a mixed order) + // Copy the graph int i = 0; for (NodeIt n(_graph); n != INVALID; ++n, ++i) { _node_id[n] = i; } - int k = std::max(int(std::sqrt(double(_arc_num))), 10); - i = 0; - for (ArcIt a(_graph); a != INVALID; ++a) { - _arc_id[a] = i; - _source[i] = _node_id[_graph.source(a)]; - _target[i] = _node_id[_graph.target(a)]; - if ((i += k) >= _arc_num) i = (i % k) + 1; + if (arc_mixing) { + // Store the arcs in a mixed order + int k = std::max(int(std::sqrt(double(_arc_num))), 10); + int i = 0, j = 0; + for (ArcIt a(_graph); a != INVALID; ++a) { + _arc_id[a] = i; + _source[i] = _node_id[_graph.source(a)]; + _target[i] = _node_id[_graph.target(a)]; + if ((i += k) >= _arc_num) i = ++j; + } + } else { + // Store the arcs in the original order + int i = 0; + for (ArcIt a(_graph); a != INVALID; ++a, ++i) { + _arc_id[a] = i; + _source[i] = _node_id[_graph.source(a)]; + _target[i] = _node_id[_graph.target(a)]; + } } - // Initialize maps - for (int i = 0; i != _node_num; ++i) { - _supply[i] = 0; - } - for (int i = 0; i != _arc_num; ++i) { - _lower[i] = 0; - _upper[i] = INF; - _cost[i] = 1; - } - _have_lower = false; - _stype = GEQ; + // Reset parameters + reset(); } /// \name Parameters @@ -768,7 +766,6 @@ /// This function sets the supply values of the nodes. /// If neither this function nor \ref stSupply() is used before /// calling \ref run(), the supply of each node will be set to zero. - /// (It makes sense only if non-zero lower bounds are given.) /// /// \param map A node map storing the supply values. /// Its \c Value type must be convertible to the \c Value type @@ -789,7 +786,6 @@ /// and the required flow value. /// If neither this function nor \ref supplyMap() is used before /// calling \ref run(), the supply of each node will be set to zero. - /// (It makes sense only if non-zero lower bounds are given.) /// /// Using this function has the same effect as using \ref supplyMap() /// with such a map in which \c k is assigned to \c s, \c -k is @@ -816,7 +812,7 @@ /// If it is not used before calling \ref run(), the \ref GEQ supply /// type will be used. /// - /// For more information see \ref SupplyType. + /// For more information, see \ref SupplyType. /// /// \return (*this) NetworkSimplex& supplyType(SupplyType supply_type) { @@ -848,11 +844,11 @@ /// that have been given are kept for the next call, unless /// \ref reset() is called, thus only the modified parameters /// have to be set again. See \ref reset() for examples. - /// However the underlying digraph must not be modified after this + /// However, the underlying digraph must not be modified after this /// class have been constructed, since it copies and extends the graph. /// /// \param pivot_rule The pivot rule that will be used during the - /// algorithm. For more information see \ref PivotRule. + /// algorithm. For more information, see \ref PivotRule. /// /// \return \c INFEASIBLE if no feasible flow exists, /// \n \c OPTIMAL if the problem has optimal solution @@ -877,7 +873,7 @@ /// It is useful for multiple run() calls. If this function is not /// used, all the parameters given before are kept for the next /// \ref run() call. - /// However the underlying digraph must not be modified after this + /// However, the underlying digraph must not be modified after this /// class have been constructed, since it copies and extends the graph. /// /// For example,