diff --git a/lemon/network_simplex.h b/lemon/network_simplex.h --- a/lemon/network_simplex.h +++ b/lemon/network_simplex.h @@ -164,9 +164,10 @@ TEMPLATE_DIGRAPH_TYPEDEFS(GR); typedef std::vector IntVector; - typedef std::vector CharVector; 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 { @@ -212,8 +213,8 @@ IntVector _succ_num; IntVector _last_succ; IntVector _dirty_revs; - CharVector _forward; - CharVector _state; + BoolVector _forward; + BoolVector _state; int _root; // Temporary data used in the current pivot iteration @@ -244,7 +245,7 @@ const IntVector &_source; const IntVector &_target; const CostVector &_cost; - const CharVector &_state; + const BoolVector &_state; const CostVector &_pi; int &_in_arc; int _search_arc_num; @@ -265,7 +266,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; @@ -273,7 +274,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; @@ -296,7 +297,7 @@ const IntVector &_source; const IntVector &_target; const CostVector &_cost; - const CharVector &_state; + const BoolVector &_state; const CostVector &_pi; int &_in_arc; int _search_arc_num; @@ -313,7 +314,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; @@ -335,7 +336,7 @@ const IntVector &_source; const IntVector &_target; const CostVector &_cost; - const CharVector &_state; + const BoolVector &_state; const CostVector &_pi; int &_in_arc; int _search_arc_num; @@ -354,7 +355,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 * @@ -367,7 +368,7 @@ Cost c, min = 0; int cnt = _block_size; int e; - 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 < min) { min = c; @@ -378,7 +379,7 @@ cnt = _block_size; } } - for (e = 0; e < _next_arc; ++e) { + for (e = 0; e != _next_arc; ++e) { c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < min) { min = c; @@ -408,7 +409,7 @@ const IntVector &_source; const IntVector &_target; const CostVector &_cost; - const CharVector &_state; + const BoolVector &_state; const CostVector &_pi; int &_in_arc; int _search_arc_num; @@ -469,7 +470,7 @@ // 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; @@ -480,7 +481,7 @@ if (_curr_length == _list_length) goto search_end; } } - for (e = 0; e < _next_arc; ++e) { + 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; @@ -511,7 +512,7 @@ const IntVector &_source; const IntVector &_target; const CostVector &_cost; - const CharVector &_state; + const BoolVector &_state; const CostVector &_pi; int &_in_arc; int _search_arc_num; @@ -564,7 +565,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]]); @@ -577,7 +578,7 @@ int cnt = _block_size; int limit = _head_length; - for (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) { @@ -589,7 +590,7 @@ cnt = _block_size; } } - for (e = 0; e < _next_arc; ++e) { + 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) { @@ -1328,7 +1329,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; } @@ -1400,6 +1401,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 @@ -1422,6 +1517,9 @@ ProblemType start() { PivotRuleImpl pivot(*this); + // Perform heuristic initial pivots + if (!initialPivots()) return UNBOUNDED; + // Execute the Network Simplex algorithm while (pivot.findEnteringArc()) { findJoinNode();