0
2
0
... | ... |
@@ -132,9 +132,9 @@ |
132 | 132 |
|
133 | 133 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
134 | 134 |
|
135 | 135 |
typedef std::vector<int> IntVector; |
136 |
typedef std::vector< |
|
136 |
typedef std::vector<char> BoolVector; |
|
137 | 137 |
typedef std::vector<Value> ValueVector; |
138 | 138 |
typedef std::vector<Cost> CostVector; |
139 | 139 |
|
140 | 140 |
private: |
... | ... |
@@ -195,8 +195,9 @@ |
195 | 195 |
{ |
196 | 196 |
private: |
197 | 197 |
|
198 | 198 |
int _node_num; |
199 |
bool _geq; |
|
199 | 200 |
const IntVector &_first_out; |
200 | 201 |
const IntVector &_target; |
201 | 202 |
const CostVector &_cost; |
202 | 203 |
const ValueVector &_res_cap; |
... | ... |
@@ -209,12 +210,12 @@ |
209 | 210 |
|
210 | 211 |
public: |
211 | 212 |
|
212 | 213 |
ResidualDijkstra(CapacityScaling& cs) : |
213 |
_node_num(cs._node_num), _first_out(cs._first_out), |
|
214 |
_target(cs._target), _cost(cs._cost), _res_cap(cs._res_cap), |
|
215 |
_excess(cs._excess), _pi(cs._pi), _pred(cs._pred), |
|
216 |
_dist(cs._node_num) |
|
214 |
_node_num(cs._node_num), _geq(cs._sum_supply < 0), |
|
215 |
_first_out(cs._first_out), _target(cs._target), _cost(cs._cost), |
|
216 |
_res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi), |
|
217 |
_pred(cs._pred), _dist(cs._node_num) |
|
217 | 218 |
{} |
218 | 219 |
|
219 | 220 |
int run(int s, Value delta = 1) { |
220 | 221 |
RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP); |
... | ... |
@@ -231,9 +232,10 @@ |
231 | 232 |
_proc_nodes.push_back(u); |
232 | 233 |
heap.pop(); |
233 | 234 |
|
234 | 235 |
// Traverse outgoing residual arcs |
235 |
|
|
236 |
int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1; |
|
237 |
for (int a = _first_out[u]; a != last_out; ++a) { |
|
236 | 238 |
if (_res_cap[a] < delta) continue; |
237 | 239 |
v = _target[a]; |
238 | 240 |
switch (heap.state(v)) { |
239 | 241 |
case Heap::PRE_HEAP: |
... | ... |
@@ -686,24 +688,27 @@ |
686 | 688 |
_sum_supply += _supply[i]; |
687 | 689 |
} |
688 | 690 |
if (_sum_supply > 0) return INFEASIBLE; |
689 | 691 |
|
690 |
// Initialize |
|
692 |
// Initialize vectors |
|
691 | 693 |
for (int i = 0; i != _root; ++i) { |
692 | 694 |
_pi[i] = 0; |
693 | 695 |
_excess[i] = _supply[i]; |
694 | 696 |
} |
695 | 697 |
|
696 | 698 |
// Remove non-zero lower bounds |
699 |
const Value MAX = std::numeric_limits<Value>::max(); |
|
700 |
int last_out; |
|
697 | 701 |
if (_have_lower) { |
698 | 702 |
for (int i = 0; i != _root; ++i) { |
699 |
|
|
703 |
last_out = _first_out[i+1]; |
|
704 |
for (int j = _first_out[i]; j != last_out; ++j) { |
|
700 | 705 |
if (_forward[j]) { |
701 | 706 |
Value c = _lower[j]; |
702 | 707 |
if (c >= 0) { |
703 |
_res_cap[j] = _upper[j] < |
|
708 |
_res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF; |
|
704 | 709 |
} else { |
705 |
_res_cap[j] = _upper[j] < |
|
710 |
_res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF; |
|
706 | 711 |
} |
707 | 712 |
_excess[i] -= c; |
708 | 713 |
_excess[_target[j]] += c; |
709 | 714 |
} else { |
... | ... |
@@ -717,17 +722,18 @@ |
717 | 722 |
} |
718 | 723 |
} |
719 | 724 |
|
720 | 725 |
// Handle negative costs |
721 |
for (int u = 0; u != _root; ++u) { |
|
722 |
for (int a = _first_out[u]; a != _first_out[u+1]; ++a) { |
|
723 |
Value rc = _res_cap[a]; |
|
724 |
if (_cost[a] < 0 && rc > 0) { |
|
725 |
if (rc == INF) return UNBOUNDED; |
|
726 |
_excess[u] -= rc; |
|
727 |
_excess[_target[a]] += rc; |
|
728 |
_res_cap[a] = 0; |
|
729 |
|
|
726 |
for (int i = 0; i != _root; ++i) { |
|
727 |
last_out = _first_out[i+1] - 1; |
|
728 |
for (int j = _first_out[i]; j != last_out; ++j) { |
|
729 |
Value rc = _res_cap[j]; |
|
730 |
if (_cost[j] < 0 && rc > 0) { |
|
731 |
if (rc >= MAX) return UNBOUNDED; |
|
732 |
_excess[i] -= rc; |
|
733 |
_excess[_target[j]] += rc; |
|
734 |
_res_cap[j] = 0; |
|
735 |
_res_cap[_reverse[j]] += rc; |
|
730 | 736 |
} |
731 | 737 |
} |
732 | 738 |
} |
733 | 739 |
|
... | ... |
@@ -735,36 +741,34 @@ |
735 | 741 |
if (_sum_supply < 0) { |
736 | 742 |
_pi[_root] = 0; |
737 | 743 |
_excess[_root] = -_sum_supply; |
738 | 744 |
for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
739 |
int u = _target[a]; |
|
740 |
if (_excess[u] < 0) { |
|
741 |
_res_cap[a] = -_excess[u] + 1; |
|
742 |
} else { |
|
743 |
_res_cap[a] = 1; |
|
744 |
} |
|
745 |
|
|
745 |
int ra = _reverse[a]; |
|
746 |
_res_cap[a] = -_sum_supply + 1; |
|
747 |
_res_cap[ra] = 0; |
|
746 | 748 |
_cost[a] = 0; |
747 |
_cost[ |
|
749 |
_cost[ra] = 0; |
|
748 | 750 |
} |
749 | 751 |
} else { |
750 | 752 |
_pi[_root] = 0; |
751 | 753 |
_excess[_root] = 0; |
752 | 754 |
for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
755 |
int ra = _reverse[a]; |
|
753 | 756 |
_res_cap[a] = 1; |
754 |
_res_cap[ |
|
757 |
_res_cap[ra] = 0; |
|
755 | 758 |
_cost[a] = 0; |
756 |
_cost[ |
|
759 |
_cost[ra] = 0; |
|
757 | 760 |
} |
758 | 761 |
} |
759 | 762 |
|
760 | 763 |
// Initialize delta value |
761 | 764 |
if (_factor > 1) { |
762 | 765 |
// With scaling |
763 | 766 |
Value max_sup = 0, max_dem = 0; |
764 | 767 |
for (int i = 0; i != _node_num; ++i) { |
765 |
if ( _excess[i] > max_sup) max_sup = _excess[i]; |
|
766 |
if (-_excess[i] > max_dem) max_dem = -_excess[i]; |
|
768 |
Value ex = _excess[i]; |
|
769 |
if ( ex > max_sup) max_sup = ex; |
|
770 |
if (-ex > max_dem) max_dem = -ex; |
|
767 | 771 |
} |
768 | 772 |
Value max_cap = 0; |
769 | 773 |
for (int j = 0; j != _res_arc_num; ++j) { |
770 | 774 |
if (_res_cap[j] > max_cap) max_cap = _res_cap[j]; |
... | ... |
@@ -788,9 +792,10 @@ |
788 | 792 |
pt = startWithoutScaling(); |
789 | 793 |
|
790 | 794 |
// Handle non-zero lower bounds |
791 | 795 |
if (_have_lower) { |
792 |
|
|
796 |
int limit = _first_out[_root]; |
|
797 |
for (int j = 0; j != limit; ++j) { |
|
793 | 798 |
if (!_forward[j]) _res_cap[j] += _lower[j]; |
794 | 799 |
} |
795 | 800 |
} |
796 | 801 |
|
... | ... |
@@ -811,10 +816,13 @@ |
811 | 816 |
int s, t; |
812 | 817 |
ResidualDijkstra _dijkstra(*this); |
813 | 818 |
while (true) { |
814 | 819 |
// Saturate all arcs not satisfying the optimality condition |
820 |
int last_out; |
|
815 | 821 |
for (int u = 0; u != _node_num; ++u) { |
816 |
|
|
822 |
last_out = _sum_supply < 0 ? |
|
823 |
_first_out[u+1] : _first_out[u+1] - 1; |
|
824 |
for (int a = _first_out[u]; a != last_out; ++a) { |
|
817 | 825 |
int v = _target[a]; |
818 | 826 |
Cost c = _cost[a] + _pi[u] - _pi[v]; |
819 | 827 |
Value rc = _res_cap[a]; |
820 | 828 |
if (c < 0 && rc >= _delta) { |
... | ... |
@@ -829,10 +837,11 @@ |
829 | 837 |
// Find excess nodes and deficit nodes |
830 | 838 |
_excess_nodes.clear(); |
831 | 839 |
_deficit_nodes.clear(); |
832 | 840 |
for (int u = 0; u != _node_num; ++u) { |
833 |
if (_excess[u] >= _delta) _excess_nodes.push_back(u); |
|
834 |
if (_excess[u] <= -_delta) _deficit_nodes.push_back(u); |
|
841 |
Value ex = _excess[u]; |
|
842 |
if (ex >= _delta) _excess_nodes.push_back(u); |
|
843 |
if (ex <= -_delta) _deficit_nodes.push_back(u); |
|
835 | 844 |
} |
836 | 845 |
int next_node = 0, next_def_node = 0; |
837 | 846 |
|
838 | 847 |
// Find augmenting shortest paths |
... | ... |
@@ -163,9 +163,9 @@ |
163 | 163 |
|
164 | 164 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
165 | 165 |
|
166 | 166 |
typedef std::vector<int> IntVector; |
167 |
typedef std::vector< |
|
167 |
typedef std::vector<char> CharVector; |
|
168 | 168 |
typedef std::vector<Value> ValueVector; |
169 | 169 |
typedef std::vector<Cost> CostVector; |
170 | 170 |
|
171 | 171 |
// State constants for arcs |
... | ... |
@@ -211,17 +211,19 @@ |
211 | 211 |
IntVector _rev_thread; |
212 | 212 |
IntVector _succ_num; |
213 | 213 |
IntVector _last_succ; |
214 | 214 |
IntVector _dirty_revs; |
215 |
BoolVector _forward; |
|
216 |
IntVector _state; |
|
215 |
CharVector _forward; |
|
216 |
CharVector _state; |
|
217 | 217 |
int _root; |
218 | 218 |
|
219 | 219 |
// Temporary data used in the current pivot iteration |
220 | 220 |
int in_arc, join, u_in, v_in, u_out, v_out; |
221 | 221 |
int first, second, right, last; |
222 | 222 |
int stem, par_stem, new_stem; |
223 | 223 |
Value delta; |
224 |
|
|
225 |
const Value MAX; |
|
224 | 226 |
|
225 | 227 |
public: |
226 | 228 |
|
227 | 229 |
/// \brief Constant for infinite upper bounds (capacities). |
... | ... |
@@ -241,9 +243,9 @@ |
241 | 243 |
// References to the NetworkSimplex class |
242 | 244 |
const IntVector &_source; |
243 | 245 |
const IntVector &_target; |
244 | 246 |
const CostVector &_cost; |
245 |
const |
|
247 |
const CharVector &_state; |
|
246 | 248 |
const CostVector &_pi; |
247 | 249 |
int &_in_arc; |
248 | 250 |
int _search_arc_num; |
249 | 251 |
|
... | ... |
@@ -293,9 +295,9 @@ |
293 | 295 |
// References to the NetworkSimplex class |
294 | 296 |
const IntVector &_source; |
295 | 297 |
const IntVector &_target; |
296 | 298 |
const CostVector &_cost; |
297 |
const |
|
299 |
const CharVector &_state; |
|
298 | 300 |
const CostVector &_pi; |
299 | 301 |
int &_in_arc; |
300 | 302 |
int _search_arc_num; |
301 | 303 |
|
... | ... |
@@ -332,9 +334,9 @@ |
332 | 334 |
// References to the NetworkSimplex class |
333 | 335 |
const IntVector &_source; |
334 | 336 |
const IntVector &_target; |
335 | 337 |
const CostVector &_cost; |
336 |
const |
|
338 |
const CharVector &_state; |
|
337 | 339 |
const CostVector &_pi; |
338 | 340 |
int &_in_arc; |
339 | 341 |
int _search_arc_num; |
340 | 342 |
|
... | ... |
@@ -405,9 +407,9 @@ |
405 | 407 |
// References to the NetworkSimplex class |
406 | 408 |
const IntVector &_source; |
407 | 409 |
const IntVector &_target; |
408 | 410 |
const CostVector &_cost; |
409 |
const |
|
411 |
const CharVector &_state; |
|
410 | 412 |
const CostVector &_pi; |
411 | 413 |
int &_in_arc; |
412 | 414 |
int _search_arc_num; |
413 | 415 |
|
... | ... |
@@ -508,9 +510,9 @@ |
508 | 510 |
// References to the NetworkSimplex class |
509 | 511 |
const IntVector &_source; |
510 | 512 |
const IntVector &_target; |
511 | 513 |
const CostVector &_cost; |
512 |
const |
|
514 |
const CharVector &_state; |
|
513 | 515 |
const CostVector &_pi; |
514 | 516 |
int &_in_arc; |
515 | 517 |
int _search_arc_num; |
516 | 518 |
|
... | ... |
@@ -630,11 +632,11 @@ |
630 | 632 |
/// In special cases, it could lead to better overall performance, |
631 | 633 |
/// but it is usually slower. Therefore it is disabled by default. |
632 | 634 |
NetworkSimplex(const GR& graph, bool arc_mixing = false) : |
633 | 635 |
_graph(graph), _node_id(graph), _arc_id(graph), |
636 |
MAX(std::numeric_limits<Value>::max()), |
|
634 | 637 |
INF(std::numeric_limits<Value>::has_infinity ? |
635 |
std::numeric_limits<Value>::infinity() : |
|
636 |
std::numeric_limits<Value>::max()) |
|
638 |
std::numeric_limits<Value>::infinity() : MAX) |
|
637 | 639 |
{ |
638 | 640 |
// Check the value types |
639 | 641 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
640 | 642 |
"The flow type of NetworkSimplex must be signed"); |
... | ... |
@@ -1019,11 +1021,11 @@ |
1019 | 1021 |
if (_have_lower) { |
1020 | 1022 |
for (int i = 0; i != _arc_num; ++i) { |
1021 | 1023 |
Value c = _lower[i]; |
1022 | 1024 |
if (c >= 0) { |
1023 |
_cap[i] = _upper[i] < |
|
1025 |
_cap[i] = _upper[i] < MAX ? _upper[i] - c : INF; |
|
1024 | 1026 |
} else { |
1025 |
_cap[i] = _upper[i] < |
|
1027 |
_cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF; |
|
1026 | 1028 |
} |
1027 | 1029 |
_supply[_source[i]] -= c; |
1028 | 1030 |
_supply[_target[i]] += c; |
1029 | 1031 |
} |
... | ... |
@@ -1213,9 +1215,9 @@ |
1213 | 1215 |
// Search the cycle along the path form the first node to the root |
1214 | 1216 |
for (int u = first; u != join; u = _parent[u]) { |
1215 | 1217 |
e = _pred[u]; |
1216 | 1218 |
d = _forward[u] ? |
1217 |
_flow[e] : (_cap[e] |
|
1219 |
_flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]); |
|
1218 | 1220 |
if (d < delta) { |
1219 | 1221 |
delta = d; |
1220 | 1222 |
u_out = u; |
1221 | 1223 |
result = 1; |
... | ... |
@@ -1224,9 +1226,9 @@ |
1224 | 1226 |
// Search the cycle along the path form the second node to the root |
1225 | 1227 |
for (int u = second; u != join; u = _parent[u]) { |
1226 | 1228 |
e = _pred[u]; |
1227 | 1229 |
d = _forward[u] ? |
1228 |
(_cap[e] |
|
1230 |
(_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e]; |
|
1229 | 1231 |
if (d <= delta) { |
1230 | 1232 |
delta = d; |
1231 | 1233 |
u_out = u; |
1232 | 1234 |
result = 2; |
... | ... |
@@ -1423,9 +1425,9 @@ |
1423 | 1425 |
// Execute the Network Simplex algorithm |
1424 | 1426 |
while (pivot.findEnteringArc()) { |
1425 | 1427 |
findJoinNode(); |
1426 | 1428 |
bool change = findLeavingArc(); |
1427 |
if (delta >= |
|
1429 |
if (delta >= MAX) return UNBOUNDED; |
|
1428 | 1430 |
changeFlow(change); |
1429 | 1431 |
if (change) { |
1430 | 1432 |
updateTreeStructure(); |
1431 | 1433 |
updatePotential(); |
0 comments (0 inline)