1.1 --- a/lemon/network_simplex.h Tue Feb 09 23:29:51 2010 +0100
1.2 +++ b/lemon/network_simplex.h Sat Feb 20 18:39:03 2010 +0100
1.3 @@ -164,9 +164,10 @@
1.4 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.5
1.6 typedef std::vector<int> IntVector;
1.7 - typedef std::vector<char> CharVector;
1.8 typedef std::vector<Value> ValueVector;
1.9 typedef std::vector<Cost> CostVector;
1.10 + typedef std::vector<char> BoolVector;
1.11 + // Note: vector<char> is used instead of vector<bool> for efficiency reasons
1.12
1.13 // State constants for arcs
1.14 enum ArcStateEnum {
1.15 @@ -212,8 +213,8 @@
1.16 IntVector _succ_num;
1.17 IntVector _last_succ;
1.18 IntVector _dirty_revs;
1.19 - CharVector _forward;
1.20 - CharVector _state;
1.21 + BoolVector _forward;
1.22 + BoolVector _state;
1.23 int _root;
1.24
1.25 // Temporary data used in the current pivot iteration
1.26 @@ -244,7 +245,7 @@
1.27 const IntVector &_source;
1.28 const IntVector &_target;
1.29 const CostVector &_cost;
1.30 - const CharVector &_state;
1.31 + const BoolVector &_state;
1.32 const CostVector &_pi;
1.33 int &_in_arc;
1.34 int _search_arc_num;
1.35 @@ -265,7 +266,7 @@
1.36 // Find next entering arc
1.37 bool findEnteringArc() {
1.38 Cost c;
1.39 - for (int e = _next_arc; e < _search_arc_num; ++e) {
1.40 + for (int e = _next_arc; e != _search_arc_num; ++e) {
1.41 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.42 if (c < 0) {
1.43 _in_arc = e;
1.44 @@ -273,7 +274,7 @@
1.45 return true;
1.46 }
1.47 }
1.48 - for (int e = 0; e < _next_arc; ++e) {
1.49 + for (int e = 0; e != _next_arc; ++e) {
1.50 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.51 if (c < 0) {
1.52 _in_arc = e;
1.53 @@ -296,7 +297,7 @@
1.54 const IntVector &_source;
1.55 const IntVector &_target;
1.56 const CostVector &_cost;
1.57 - const CharVector &_state;
1.58 + const BoolVector &_state;
1.59 const CostVector &_pi;
1.60 int &_in_arc;
1.61 int _search_arc_num;
1.62 @@ -313,7 +314,7 @@
1.63 // Find next entering arc
1.64 bool findEnteringArc() {
1.65 Cost c, min = 0;
1.66 - for (int e = 0; e < _search_arc_num; ++e) {
1.67 + for (int e = 0; e != _search_arc_num; ++e) {
1.68 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.69 if (c < min) {
1.70 min = c;
1.71 @@ -335,7 +336,7 @@
1.72 const IntVector &_source;
1.73 const IntVector &_target;
1.74 const CostVector &_cost;
1.75 - const CharVector &_state;
1.76 + const BoolVector &_state;
1.77 const CostVector &_pi;
1.78 int &_in_arc;
1.79 int _search_arc_num;
1.80 @@ -354,7 +355,7 @@
1.81 _next_arc(0)
1.82 {
1.83 // The main parameters of the pivot rule
1.84 - const double BLOCK_SIZE_FACTOR = 0.5;
1.85 + const double BLOCK_SIZE_FACTOR = 1.0;
1.86 const int MIN_BLOCK_SIZE = 10;
1.87
1.88 _block_size = std::max( int(BLOCK_SIZE_FACTOR *
1.89 @@ -367,7 +368,7 @@
1.90 Cost c, min = 0;
1.91 int cnt = _block_size;
1.92 int e;
1.93 - for (e = _next_arc; e < _search_arc_num; ++e) {
1.94 + for (e = _next_arc; e != _search_arc_num; ++e) {
1.95 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.96 if (c < min) {
1.97 min = c;
1.98 @@ -378,7 +379,7 @@
1.99 cnt = _block_size;
1.100 }
1.101 }
1.102 - for (e = 0; e < _next_arc; ++e) {
1.103 + for (e = 0; e != _next_arc; ++e) {
1.104 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.105 if (c < min) {
1.106 min = c;
1.107 @@ -408,7 +409,7 @@
1.108 const IntVector &_source;
1.109 const IntVector &_target;
1.110 const CostVector &_cost;
1.111 - const CharVector &_state;
1.112 + const BoolVector &_state;
1.113 const CostVector &_pi;
1.114 int &_in_arc;
1.115 int _search_arc_num;
1.116 @@ -469,7 +470,7 @@
1.117 // Major iteration: build a new candidate list
1.118 min = 0;
1.119 _curr_length = 0;
1.120 - for (e = _next_arc; e < _search_arc_num; ++e) {
1.121 + for (e = _next_arc; e != _search_arc_num; ++e) {
1.122 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.123 if (c < 0) {
1.124 _candidates[_curr_length++] = e;
1.125 @@ -480,7 +481,7 @@
1.126 if (_curr_length == _list_length) goto search_end;
1.127 }
1.128 }
1.129 - for (e = 0; e < _next_arc; ++e) {
1.130 + for (e = 0; e != _next_arc; ++e) {
1.131 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.132 if (c < 0) {
1.133 _candidates[_curr_length++] = e;
1.134 @@ -511,7 +512,7 @@
1.135 const IntVector &_source;
1.136 const IntVector &_target;
1.137 const CostVector &_cost;
1.138 - const CharVector &_state;
1.139 + const BoolVector &_state;
1.140 const CostVector &_pi;
1.141 int &_in_arc;
1.142 int _search_arc_num;
1.143 @@ -564,7 +565,7 @@
1.144 bool findEnteringArc() {
1.145 // Check the current candidate list
1.146 int e;
1.147 - for (int i = 0; i < _curr_length; ++i) {
1.148 + for (int i = 0; i != _curr_length; ++i) {
1.149 e = _candidates[i];
1.150 _cand_cost[e] = _state[e] *
1.151 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.152 @@ -577,7 +578,7 @@
1.153 int cnt = _block_size;
1.154 int limit = _head_length;
1.155
1.156 - for (e = _next_arc; e < _search_arc_num; ++e) {
1.157 + for (e = _next_arc; e != _search_arc_num; ++e) {
1.158 _cand_cost[e] = _state[e] *
1.159 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.160 if (_cand_cost[e] < 0) {
1.161 @@ -589,7 +590,7 @@
1.162 cnt = _block_size;
1.163 }
1.164 }
1.165 - for (e = 0; e < _next_arc; ++e) {
1.166 + for (e = 0; e != _next_arc; ++e) {
1.167 _cand_cost[e] = _state[e] *
1.168 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.169 if (_cand_cost[e] < 0) {
1.170 @@ -1328,7 +1329,7 @@
1.171 }
1.172
1.173 // Update _rev_thread using the new _thread values
1.174 - for (int i = 0; i < int(_dirty_revs.size()); ++i) {
1.175 + for (int i = 0; i != int(_dirty_revs.size()); ++i) {
1.176 u = _dirty_revs[i];
1.177 _rev_thread[_thread[u]] = u;
1.178 }
1.179 @@ -1400,6 +1401,100 @@
1.180 }
1.181 }
1.182
1.183 + // Heuristic initial pivots
1.184 + bool initialPivots() {
1.185 + Value curr, total = 0;
1.186 + std::vector<Node> supply_nodes, demand_nodes;
1.187 + for (NodeIt u(_graph); u != INVALID; ++u) {
1.188 + curr = _supply[_node_id[u]];
1.189 + if (curr > 0) {
1.190 + total += curr;
1.191 + supply_nodes.push_back(u);
1.192 + }
1.193 + else if (curr < 0) {
1.194 + demand_nodes.push_back(u);
1.195 + }
1.196 + }
1.197 + if (_sum_supply > 0) total -= _sum_supply;
1.198 + if (total <= 0) return true;
1.199 +
1.200 + IntVector arc_vector;
1.201 + if (_sum_supply >= 0) {
1.202 + if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {
1.203 + // Perform a reverse graph search from the sink to the source
1.204 + typename GR::template NodeMap<bool> reached(_graph, false);
1.205 + Node s = supply_nodes[0], t = demand_nodes[0];
1.206 + std::vector<Node> stack;
1.207 + reached[t] = true;
1.208 + stack.push_back(t);
1.209 + while (!stack.empty()) {
1.210 + Node u, v = stack.back();
1.211 + stack.pop_back();
1.212 + if (v == s) break;
1.213 + for (InArcIt a(_graph, v); a != INVALID; ++a) {
1.214 + if (reached[u = _graph.source(a)]) continue;
1.215 + int j = _arc_id[a];
1.216 + if (_cap[j] >= total) {
1.217 + arc_vector.push_back(j);
1.218 + reached[u] = true;
1.219 + stack.push_back(u);
1.220 + }
1.221 + }
1.222 + }
1.223 + } else {
1.224 + // Find the min. cost incomming arc for each demand node
1.225 + for (int i = 0; i != int(demand_nodes.size()); ++i) {
1.226 + Node v = demand_nodes[i];
1.227 + Cost c, min_cost = std::numeric_limits<Cost>::max();
1.228 + Arc min_arc = INVALID;
1.229 + for (InArcIt a(_graph, v); a != INVALID; ++a) {
1.230 + c = _cost[_arc_id[a]];
1.231 + if (c < min_cost) {
1.232 + min_cost = c;
1.233 + min_arc = a;
1.234 + }
1.235 + }
1.236 + if (min_arc != INVALID) {
1.237 + arc_vector.push_back(_arc_id[min_arc]);
1.238 + }
1.239 + }
1.240 + }
1.241 + } else {
1.242 + // Find the min. cost outgoing arc for each supply node
1.243 + for (int i = 0; i != int(supply_nodes.size()); ++i) {
1.244 + Node u = supply_nodes[i];
1.245 + Cost c, min_cost = std::numeric_limits<Cost>::max();
1.246 + Arc min_arc = INVALID;
1.247 + for (OutArcIt a(_graph, u); a != INVALID; ++a) {
1.248 + c = _cost[_arc_id[a]];
1.249 + if (c < min_cost) {
1.250 + min_cost = c;
1.251 + min_arc = a;
1.252 + }
1.253 + }
1.254 + if (min_arc != INVALID) {
1.255 + arc_vector.push_back(_arc_id[min_arc]);
1.256 + }
1.257 + }
1.258 + }
1.259 +
1.260 + // Perform heuristic initial pivots
1.261 + for (int i = 0; i != int(arc_vector.size()); ++i) {
1.262 + in_arc = arc_vector[i];
1.263 + if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -
1.264 + _pi[_target[in_arc]]) >= 0) continue;
1.265 + findJoinNode();
1.266 + bool change = findLeavingArc();
1.267 + if (delta >= MAX) return false;
1.268 + changeFlow(change);
1.269 + if (change) {
1.270 + updateTreeStructure();
1.271 + updatePotential();
1.272 + }
1.273 + }
1.274 + return true;
1.275 + }
1.276 +
1.277 // Execute the algorithm
1.278 ProblemType start(PivotRule pivot_rule) {
1.279 // Select the pivot rule implementation
1.280 @@ -1422,6 +1517,9 @@
1.281 ProblemType start() {
1.282 PivotRuleImpl pivot(*this);
1.283
1.284 + // Perform heuristic initial pivots
1.285 + if (!initialPivots()) return UNBOUNDED;
1.286 +
1.287 // Execute the Network Simplex algorithm
1.288 while (pivot.findEnteringArc()) {
1.289 findJoinNode();