diff --git a/lemon/network_simplex.h b/lemon/network_simplex.h --- a/lemon/network_simplex.h +++ b/lemon/network_simplex.h @@ -2,7 +2,7 @@ * * This file is a part of LEMON, a generic C++ optimization library. * - * Copyright (C) 2003-2009 + * Copyright (C) 2003-2010 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport * (Egervary Research Group on Combinatorial Optimization, EGRES). * @@ -40,15 +40,17 @@ /// 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". - /// 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. + /// for finding a \ref min_cost_flow "minimum cost flow" + /// \ref amo93networkflows, \ref dantzig63linearprog, + /// \ref kellyoneill91netsimplex. + /// This algorithm is a highly efficient specialized version of the + /// linear programming simplex method directly for the minimum cost + /// flow problem. /// - /// 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. + /// In general, %NetworkSimplex is the fastest implementation available + /// in LEMON for this problem. + /// Moreover, it supports both directions of the supply/demand inequality + /// 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 @@ -56,17 +58,17 @@ /// specified, then default values will be used. /// /// \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. - /// \tparam C The value type used for costs and potentials in the - /// algorithm. By default it is the same as \c V. + /// \tparam V The number type used for flow amounts, capacity bounds + /// and supply values in the algorithm. By default, it is \c int. + /// \tparam C The number type used for costs and potentials in the + /// algorithm. By default, it is the same as \c V. /// - /// \warning Both value types must be signed and all input data must + /// \warning Both number 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 { @@ -95,7 +97,7 @@ /// infinite upper bound. UNBOUNDED }; - + /// \brief Constants for selecting the type of the supply constraints. /// /// Enum type containing constants for selecting the supply type, @@ -113,7 +115,7 @@ /// supply/demand constraints in the definition of the problem. LEQ }; - + /// \brief Constants for selecting the pivot rule. /// /// Enum type containing constants for selecting the pivot rule for @@ -122,59 +124,62 @@ /// \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() + /// test inputs. + /// 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. ALTERING_LIST }; - + private: TEMPLATE_DIGRAPH_TYPEDEFS(GR); - typedef std::vector ArcVector; - typedef std::vector NodeVector; typedef std::vector IntVector; - typedef std::vector BoolVector; typedef std::vector ValueVector; typedef std::vector CostVector; + typedef std::vector BoolVector; + // Note: vector is used instead of vector for efficiency reasons // State constants for arcs - enum ArcStateEnum { + enum ArcState { STATE_UPPER = -1, STATE_TREE = 0, STATE_LOWER = 1 }; + typedef std::vector StateVector; + // Note: vector is used instead of vector for + // efficiency reasons + private: // Data related to the underlying digraph @@ -194,6 +199,7 @@ IntArcMap _arc_id; IntVector _source; IntVector _target; + bool _arc_mixing; // Node and arc data ValueVector _lower; @@ -213,7 +219,7 @@ IntVector _last_succ; IntVector _dirty_revs; BoolVector _forward; - IntVector _state; + StateVector _state; int _root; // Temporary data used in the current pivot iteration @@ -222,8 +228,10 @@ int stem, par_stem, new_stem; Value delta; + const Value MAX; + public: - + /// \brief Constant for infinite upper bounds (capacities). /// /// Constant for infinite upper bounds (capacities). @@ -242,7 +250,7 @@ const IntVector &_source; const IntVector &_target; const CostVector &_cost; - const IntVector &_state; + const StateVector &_state; const CostVector &_pi; int &_in_arc; int _search_arc_num; @@ -263,7 +271,7 @@ // Find next entering arc bool findEnteringArc() { Cost c; - for (int e = _next_arc; e < _search_arc_num; ++e) { + for (int e = _next_arc; e != _search_arc_num; ++e) { c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < 0) { _in_arc = e; @@ -271,7 +279,7 @@ return true; } } - for (int e = 0; e < _next_arc; ++e) { + for (int e = 0; e != _next_arc; ++e) { c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < 0) { _in_arc = e; @@ -294,7 +302,7 @@ const IntVector &_source; const IntVector &_target; const CostVector &_cost; - const IntVector &_state; + const StateVector &_state; const CostVector &_pi; int &_in_arc; int _search_arc_num; @@ -311,7 +319,7 @@ // Find next entering arc bool findEnteringArc() { Cost c, min = 0; - for (int e = 0; e < _search_arc_num; ++e) { + for (int e = 0; e != _search_arc_num; ++e) { c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < min) { min = c; @@ -333,7 +341,7 @@ const IntVector &_source; const IntVector &_target; const CostVector &_cost; - const IntVector &_state; + const StateVector &_state; const CostVector &_pi; int &_in_arc; int _search_arc_num; @@ -352,7 +360,7 @@ _next_arc(0) { // The main parameters of the pivot rule - const double BLOCK_SIZE_FACTOR = 0.5; + const double BLOCK_SIZE_FACTOR = 1.0; const int MIN_BLOCK_SIZE = 10; _block_size = std::max( int(BLOCK_SIZE_FACTOR * @@ -364,33 +372,32 @@ bool findEnteringArc() { Cost c, min = 0; int cnt = _block_size; - int e, min_arc = _next_arc; - for (e = _next_arc; e < _search_arc_num; ++e) { + 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; } @@ -407,7 +414,7 @@ const IntVector &_source; const IntVector &_target; const CostVector &_cost; - const IntVector &_state; + const StateVector &_state; const CostVector &_pi; int &_in_arc; int _search_arc_num; @@ -428,7 +435,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 +452,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,48 +463,44 @@ 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 min = 0; _curr_length = 0; - for (e = _next_arc; e < _search_arc_num; ++e) { + for (e = _next_arc; e != _search_arc_num; ++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; + _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; } @@ -514,7 +517,7 @@ const IntVector &_source; const IntVector &_target; const CostVector &_cost; - const IntVector &_state; + const StateVector &_state; const CostVector &_pi; int &_in_arc; int _search_arc_num; @@ -549,7 +552,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; @@ -567,7 +570,7 @@ bool findEnteringArc() { // Check the current candidate list int e; - for (int i = 0; i < _curr_length; ++i) { + for (int i = 0; i != _curr_length; ++i) { e = _candidates[i]; _cand_cost[e] = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); @@ -578,39 +581,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 +617,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,69 +633,25 @@ /// 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), + _arc_mixing(arc_mixing), + MAX(std::numeric_limits::max()), INF(std::numeric_limits::has_infinity ? - std::numeric_limits::infinity() : - std::numeric_limits::max()) + std::numeric_limits::infinity() : MAX) { - // Check the value types + // Check the number types LEMON_ASSERT(std::numeric_limits::is_signed, "The flow type of NetworkSimplex must be signed"); LEMON_ASSERT(std::numeric_limits::is_signed, "The cost type of NetworkSimplex must be signed"); - - // Resize vectors - _node_num = countNodes(_graph); - _arc_num = countArcs(_graph); - int all_node_num = _node_num + 1; - int max_arc_num = _arc_num + 2 * _node_num; - _source.resize(max_arc_num); - _target.resize(max_arc_num); - - _lower.resize(_arc_num); - _upper.resize(_arc_num); - _cap.resize(max_arc_num); - _cost.resize(max_arc_num); - _supply.resize(all_node_num); - _flow.resize(max_arc_num); - _pi.resize(all_node_num); - - _parent.resize(all_node_num); - _pred.resize(all_node_num); - _forward.resize(all_node_num); - _thread.resize(all_node_num); - _rev_thread.resize(all_node_num); - _succ_num.resize(all_node_num); - _last_succ.resize(all_node_num); - _state.resize(max_arc_num); - - // Copy the graph (store the arcs in a mixed order) - 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; - } - - // 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 data structures + reset(); } /// \name Parameters @@ -729,7 +685,7 @@ /// This function sets the upper bounds (capacities) on the arcs. /// If it is not used before calling \ref run(), the upper bounds /// will be set to \ref INF on all arcs (i.e. the flow value will be - /// unbounded from above on each arc). + /// unbounded from above). /// /// \param map An arc map storing the upper bounds. /// Its \c Value type must be convertible to the \c Value type @@ -768,7 +724,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 +744,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 @@ -809,14 +763,14 @@ _supply[_node_id[t]] = -k; return *this; } - + /// \brief Set the type of the supply constraints. /// /// This function sets the type of the supply/demand constraints. /// 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) { @@ -835,7 +789,7 @@ /// /// This function runs the algorithm. /// The paramters can be specified using functions \ref lowerMap(), - /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), /// \ref supplyType(). /// For example, /// \code @@ -844,15 +798,15 @@ /// .supplyMap(sup).run(); /// \endcode /// - /// This function can be called more than once. All the parameters - /// 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 - /// class have been constructed, since it copies and extends the graph. + /// This function can be called more than once. All the given parameters + /// are kept for the next call, unless \ref resetParams() or \ref reset() + /// is used, thus only the modified parameters have to be set again. + /// If the underlying digraph was also modified after the construction + /// of the class (or the last \ref reset() call), then the \ref reset() + /// function must be called. /// /// \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 @@ -863,6 +817,7 @@ /// cost and infinite upper bound. /// /// \see ProblemType, PivotRule + /// \see resetParams(), reset() ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) { if (!init()) return INFEASIBLE; return start(pivot_rule); @@ -874,11 +829,12 @@ /// before using functions \ref lowerMap(), \ref upperMap(), /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(). /// - /// 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 - /// class have been constructed, since it copies and extends the graph. + /// It is useful for multiple \ref run() calls. Basically, all the given + /// parameters are kept for the next \ref run() call, unless + /// \ref resetParams() or \ref reset() is used. + /// If the underlying digraph was also modified after the construction + /// of the class or the last \ref reset() call, then the \ref reset() + /// function must be used, otherwise \ref resetParams() is sufficient. /// /// For example, /// \code @@ -888,20 +844,22 @@ /// ns.lowerMap(lower).upperMap(upper).costMap(cost) /// .supplyMap(sup).run(); /// - /// // Run again with modified cost map (reset() is not called, + /// // Run again with modified cost map (resetParams() is not called, /// // so only the cost map have to be set again) /// cost[e] += 100; /// ns.costMap(cost).run(); /// - /// // Run again from scratch using reset() + /// // Run again from scratch using resetParams() /// // (the lower bounds will be set to zero on all arcs) - /// ns.reset(); + /// ns.resetParams(); /// ns.upperMap(capacity).costMap(cost) /// .supplyMap(sup).run(); /// \endcode /// /// \return (*this) - NetworkSimplex& reset() { + /// + /// \see reset(), run() + NetworkSimplex& resetParams() { for (int i = 0; i != _node_num; ++i) { _supply[i] = 0; } @@ -915,6 +873,83 @@ return *this; } + /// \brief Reset the internal data structures and all the parameters + /// that have been given before. + /// + /// This function resets the internal data structures and all the + /// paramaters that have been given before using functions \ref lowerMap(), + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), + /// \ref supplyType(). + /// + /// It is useful for multiple \ref run() calls. Basically, all the given + /// parameters are kept for the next \ref run() call, unless + /// \ref resetParams() or \ref reset() is used. + /// If the underlying digraph was also modified after the construction + /// of the class or the last \ref reset() call, then the \ref reset() + /// function must be used, otherwise \ref resetParams() is sufficient. + /// + /// See \ref resetParams() for examples. + /// + /// \return (*this) + /// + /// \see resetParams(), run() + NetworkSimplex& reset() { + // Resize vectors + _node_num = countNodes(_graph); + _arc_num = countArcs(_graph); + int all_node_num = _node_num + 1; + int max_arc_num = _arc_num + 2 * _node_num; + + _source.resize(max_arc_num); + _target.resize(max_arc_num); + + _lower.resize(_arc_num); + _upper.resize(_arc_num); + _cap.resize(max_arc_num); + _cost.resize(max_arc_num); + _supply.resize(all_node_num); + _flow.resize(max_arc_num); + _pi.resize(all_node_num); + + _parent.resize(all_node_num); + _pred.resize(all_node_num); + _forward.resize(all_node_num); + _thread.resize(all_node_num); + _rev_thread.resize(all_node_num); + _succ_num.resize(all_node_num); + _last_succ.resize(all_node_num); + _state.resize(max_arc_num); + + // Copy the graph + int i = 0; + for (NodeIt n(_graph); n != INVALID; ++n, ++i) { + _node_id[n] = i; + } + 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)]; + } + } + + // Reset parameters + resetParams(); + return *this; + } + /// @} /// \name Query Functions @@ -1024,9 +1059,9 @@ for (int i = 0; i != _arc_num; ++i) { Value c = _lower[i]; if (c >= 0) { - _cap[i] = _upper[i] < INF ? _upper[i] - c : INF; + _cap[i] = _upper[i] < MAX ? _upper[i] - c : INF; } else { - _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF; + _cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF; } _supply[_source[i]] -= c; _supply[_target[i]] += c; @@ -1054,7 +1089,7 @@ _flow[i] = 0; _state[i] = STATE_LOWER; } - + // Set data for the artificial root node _root = _node_num; _parent[_root] = -1; @@ -1218,7 +1253,7 @@ for (int u = first; u != join; u = _parent[u]) { e = _pred[u]; d = _forward[u] ? - _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]); + _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]); if (d < delta) { delta = d; u_out = u; @@ -1228,8 +1263,8 @@ // Search the cycle along the path form the second node to the root for (int u = second; u != join; u = _parent[u]) { e = _pred[u]; - d = _forward[u] ? - (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e]; + d = _forward[u] ? + (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e]; if (d <= delta) { delta = d; u_out = u; @@ -1330,7 +1365,7 @@ } // Update _rev_thread using the new _thread values - for (int i = 0; i < int(_dirty_revs.size()); ++i) { + for (int i = 0; i != int(_dirty_revs.size()); ++i) { u = _dirty_revs[i]; _rev_thread[_thread[u]] = u; } @@ -1402,6 +1437,100 @@ } } + // Heuristic initial pivots + bool initialPivots() { + Value curr, total = 0; + std::vector supply_nodes, demand_nodes; + for (NodeIt u(_graph); u != INVALID; ++u) { + curr = _supply[_node_id[u]]; + if (curr > 0) { + total += curr; + supply_nodes.push_back(u); + } + else if (curr < 0) { + demand_nodes.push_back(u); + } + } + if (_sum_supply > 0) total -= _sum_supply; + if (total <= 0) return true; + + IntVector arc_vector; + if (_sum_supply >= 0) { + if (supply_nodes.size() == 1 && demand_nodes.size() == 1) { + // Perform a reverse graph search from the sink to the source + typename GR::template NodeMap reached(_graph, false); + Node s = supply_nodes[0], t = demand_nodes[0]; + std::vector stack; + reached[t] = true; + stack.push_back(t); + while (!stack.empty()) { + Node u, v = stack.back(); + stack.pop_back(); + if (v == s) break; + for (InArcIt a(_graph, v); a != INVALID; ++a) { + if (reached[u = _graph.source(a)]) continue; + int j = _arc_id[a]; + if (_cap[j] >= total) { + arc_vector.push_back(j); + reached[u] = true; + stack.push_back(u); + } + } + } + } else { + // Find the min. cost incomming arc for each demand node + for (int i = 0; i != int(demand_nodes.size()); ++i) { + Node v = demand_nodes[i]; + Cost c, min_cost = std::numeric_limits::max(); + Arc min_arc = INVALID; + for (InArcIt a(_graph, v); a != INVALID; ++a) { + c = _cost[_arc_id[a]]; + if (c < min_cost) { + min_cost = c; + min_arc = a; + } + } + if (min_arc != INVALID) { + arc_vector.push_back(_arc_id[min_arc]); + } + } + } + } else { + // Find the min. cost outgoing arc for each supply node + for (int i = 0; i != int(supply_nodes.size()); ++i) { + Node u = supply_nodes[i]; + Cost c, min_cost = std::numeric_limits::max(); + Arc min_arc = INVALID; + for (OutArcIt a(_graph, u); a != INVALID; ++a) { + c = _cost[_arc_id[a]]; + if (c < min_cost) { + min_cost = c; + min_arc = a; + } + } + if (min_arc != INVALID) { + arc_vector.push_back(_arc_id[min_arc]); + } + } + } + + // Perform heuristic initial pivots + for (int i = 0; i != int(arc_vector.size()); ++i) { + in_arc = arc_vector[i]; + if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] - + _pi[_target[in_arc]]) >= 0) continue; + findJoinNode(); + bool change = findLeavingArc(); + if (delta >= MAX) return false; + changeFlow(change); + if (change) { + updateTreeStructure(); + updatePotential(); + } + } + return true; + } + // Execute the algorithm ProblemType start(PivotRule pivot_rule) { // Select the pivot rule implementation @@ -1424,18 +1553,21 @@ ProblemType start() { PivotRuleImpl pivot(*this); + // Perform heuristic initial pivots + if (!initialPivots()) return UNBOUNDED; + // Execute the Network Simplex algorithm while (pivot.findEnteringArc()) { findJoinNode(); bool change = findLeavingArc(); - if (delta >= INF) return UNBOUNDED; + if (delta >= MAX) return UNBOUNDED; changeFlow(change); if (change) { updateTreeStructure(); updatePotential(); } } - + // Check feasibility for (int e = _search_arc_num; e != _all_arc_num; ++e) { if (_flow[e] != 0) return INFEASIBLE; @@ -1452,7 +1584,7 @@ } } } - + // Shift potentials to meet the requirements of the GEQ/LEQ type // optimality conditions if (_sum_supply == 0) {