... | ... |
@@ -866,97 +866,97 @@ |
866 | 866 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
867 | 867 |
_excess[_node_id[n]] = sup[n]; |
868 | 868 |
} |
869 | 869 |
for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
870 | 870 |
int u = _target[a]; |
871 | 871 |
int ra = _reverse[a]; |
872 | 872 |
_res_cap[a] = -_sum_supply + 1; |
873 | 873 |
_res_cap[ra] = -_excess[u]; |
874 | 874 |
_cost[a] = 0; |
875 | 875 |
_cost[ra] = 0; |
876 | 876 |
_excess[u] = 0; |
877 | 877 |
} |
878 | 878 |
} else { |
879 | 879 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
880 | 880 |
Value fa = flow[a]; |
881 | 881 |
_res_cap[_arc_idf[a]] = cap[a] - fa; |
882 | 882 |
_res_cap[_arc_idb[a]] = fa; |
883 | 883 |
} |
884 | 884 |
for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
885 | 885 |
int ra = _reverse[a]; |
886 | 886 |
_res_cap[a] = 0; |
887 | 887 |
_res_cap[ra] = 0; |
888 | 888 |
_cost[a] = 0; |
889 | 889 |
_cost[ra] = 0; |
890 | 890 |
} |
891 | 891 |
} |
892 | 892 |
|
893 | 893 |
return OPTIMAL; |
894 | 894 |
} |
895 | 895 |
|
896 | 896 |
// Execute the algorithm and transform the results |
897 | 897 |
void start(Method method) { |
898 | 898 |
// Maximum path length for partial augment |
899 | 899 |
const int MAX_PATH_LENGTH = 4; |
900 | 900 |
|
901 | 901 |
// Initialize data structures for buckets |
902 | 902 |
_max_rank = _alpha * _res_node_num; |
903 | 903 |
_buckets.resize(_max_rank); |
904 | 904 |
_bucket_next.resize(_res_node_num + 1); |
905 | 905 |
_bucket_prev.resize(_res_node_num + 1); |
906 | 906 |
_rank.resize(_res_node_num + 1); |
907 | 907 |
|
908 | 908 |
// Execute the algorithm |
909 | 909 |
switch (method) { |
910 | 910 |
case PUSH: |
911 | 911 |
startPush(); |
912 | 912 |
break; |
913 | 913 |
case AUGMENT: |
914 |
startAugment(); |
|
914 |
startAugment(_res_node_num - 1); |
|
915 | 915 |
break; |
916 | 916 |
case PARTIAL_AUGMENT: |
917 | 917 |
startAugment(MAX_PATH_LENGTH); |
918 | 918 |
break; |
919 | 919 |
} |
920 | 920 |
|
921 | 921 |
// Compute node potentials for the original costs |
922 | 922 |
_arc_vec.clear(); |
923 | 923 |
_cost_vec.clear(); |
924 | 924 |
for (int j = 0; j != _res_arc_num; ++j) { |
925 | 925 |
if (_res_cap[j] > 0) { |
926 | 926 |
_arc_vec.push_back(IntPair(_source[j], _target[j])); |
927 | 927 |
_cost_vec.push_back(_scost[j]); |
928 | 928 |
} |
929 | 929 |
} |
930 | 930 |
_sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); |
931 | 931 |
|
932 | 932 |
typename BellmanFord<StaticDigraph, LargeCostArcMap> |
933 | 933 |
::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map); |
934 | 934 |
bf.distMap(_pi_map); |
935 | 935 |
bf.init(0); |
936 | 936 |
bf.start(); |
937 | 937 |
|
938 | 938 |
// Handle non-zero lower bounds |
939 | 939 |
if (_have_lower) { |
940 | 940 |
int limit = _first_out[_root]; |
941 | 941 |
for (int j = 0; j != limit; ++j) { |
942 | 942 |
if (!_forward[j]) _res_cap[j] += _lower[j]; |
943 | 943 |
} |
944 | 944 |
} |
945 | 945 |
} |
946 | 946 |
|
947 | 947 |
// Initialize a cost scaling phase |
948 | 948 |
void initPhase() { |
949 | 949 |
// Saturate arcs not satisfying the optimality condition |
950 | 950 |
for (int u = 0; u != _res_node_num; ++u) { |
951 | 951 |
int last_out = _first_out[u+1]; |
952 | 952 |
LargeCost pi_u = _pi[u]; |
953 | 953 |
for (int a = _first_out[u]; a != last_out; ++a) { |
954 | 954 |
int v = _target[a]; |
955 | 955 |
if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) { |
956 | 956 |
Value delta = _res_cap[a]; |
957 | 957 |
_excess[u] -= delta; |
958 | 958 |
_excess[v] += delta; |
959 | 959 |
_res_cap[a] = 0; |
960 | 960 |
_res_cap[_reverse[a]] += delta; |
961 | 961 |
} |
962 | 962 |
} |
... | ... |
@@ -1043,97 +1043,97 @@ |
1043 | 1043 |
int new_rank_v = old_rank_v; |
1044 | 1044 |
if (nrc < LargeCost(_max_rank)) |
1045 | 1045 |
new_rank_v = r + 1 + int(nrc); |
1046 | 1046 |
|
1047 | 1047 |
// Change the rank of v |
1048 | 1048 |
if (new_rank_v < old_rank_v) { |
1049 | 1049 |
_rank[v] = new_rank_v; |
1050 | 1050 |
_next_out[v] = _first_out[v]; |
1051 | 1051 |
|
1052 | 1052 |
// Remove v from its old bucket |
1053 | 1053 |
if (old_rank_v < _max_rank) { |
1054 | 1054 |
if (_buckets[old_rank_v] == v) { |
1055 | 1055 |
_buckets[old_rank_v] = _bucket_next[v]; |
1056 | 1056 |
} else { |
1057 | 1057 |
_bucket_next[_bucket_prev[v]] = _bucket_next[v]; |
1058 | 1058 |
_bucket_prev[_bucket_next[v]] = _bucket_prev[v]; |
1059 | 1059 |
} |
1060 | 1060 |
} |
1061 | 1061 |
|
1062 | 1062 |
// Insert v to its new bucket |
1063 | 1063 |
_bucket_next[v] = _buckets[new_rank_v]; |
1064 | 1064 |
_bucket_prev[_buckets[new_rank_v]] = v; |
1065 | 1065 |
_buckets[new_rank_v] = v; |
1066 | 1066 |
} |
1067 | 1067 |
} |
1068 | 1068 |
} |
1069 | 1069 |
} |
1070 | 1070 |
|
1071 | 1071 |
// Finish search if there are no more active nodes |
1072 | 1072 |
if (_excess[u] > 0) { |
1073 | 1073 |
total_excess -= _excess[u]; |
1074 | 1074 |
if (total_excess <= 0) break; |
1075 | 1075 |
} |
1076 | 1076 |
} |
1077 | 1077 |
if (total_excess <= 0) break; |
1078 | 1078 |
} |
1079 | 1079 |
|
1080 | 1080 |
// Relabel nodes |
1081 | 1081 |
for (int u = 0; u != _res_node_num; ++u) { |
1082 | 1082 |
int k = std::min(_rank[u], r); |
1083 | 1083 |
if (k > 0) { |
1084 | 1084 |
_pi[u] -= _epsilon * k; |
1085 | 1085 |
_next_out[u] = _first_out[u]; |
1086 | 1086 |
} |
1087 | 1087 |
} |
1088 | 1088 |
} |
1089 | 1089 |
|
1090 | 1090 |
/// Execute the algorithm performing augment and relabel operations |
1091 |
void startAugment(int max_length |
|
1091 |
void startAugment(int max_length) { |
|
1092 | 1092 |
// Paramters for heuristics |
1093 | 1093 |
const int EARLY_TERM_EPSILON_LIMIT = 1000; |
1094 | 1094 |
const double GLOBAL_UPDATE_FACTOR = 3.0; |
1095 | 1095 |
|
1096 | 1096 |
const int global_update_freq = int(GLOBAL_UPDATE_FACTOR * |
1097 | 1097 |
(_res_node_num + _sup_node_num * _sup_node_num)); |
1098 | 1098 |
int next_update_limit = global_update_freq; |
1099 | 1099 |
|
1100 | 1100 |
int relabel_cnt = 0; |
1101 | 1101 |
|
1102 | 1102 |
// Perform cost scaling phases |
1103 | 1103 |
std::vector<int> path; |
1104 | 1104 |
for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? |
1105 | 1105 |
1 : _epsilon / _alpha ) |
1106 | 1106 |
{ |
1107 | 1107 |
// Early termination heuristic |
1108 | 1108 |
if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { |
1109 | 1109 |
if (earlyTermination()) break; |
1110 | 1110 |
} |
1111 | 1111 |
|
1112 | 1112 |
// Initialize current phase |
1113 | 1113 |
initPhase(); |
1114 | 1114 |
|
1115 | 1115 |
// Perform partial augment and relabel operations |
1116 | 1116 |
while (true) { |
1117 | 1117 |
// Select an active node (FIFO selection) |
1118 | 1118 |
while (_active_nodes.size() > 0 && |
1119 | 1119 |
_excess[_active_nodes.front()] <= 0) { |
1120 | 1120 |
_active_nodes.pop_front(); |
1121 | 1121 |
} |
1122 | 1122 |
if (_active_nodes.size() == 0) break; |
1123 | 1123 |
int start = _active_nodes.front(); |
1124 | 1124 |
|
1125 | 1125 |
// Find an augmenting path from the start node |
1126 | 1126 |
path.clear(); |
1127 | 1127 |
int tip = start; |
1128 | 1128 |
while (_excess[tip] >= 0 && int(path.size()) < max_length) { |
1129 | 1129 |
int u; |
1130 | 1130 |
LargeCost min_red_cost, rc, pi_tip = _pi[tip]; |
1131 | 1131 |
int last_out = _first_out[tip+1]; |
1132 | 1132 |
for (int a = _next_out[tip]; a != last_out; ++a) { |
1133 | 1133 |
u = _target[a]; |
1134 | 1134 |
if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) { |
1135 | 1135 |
path.push_back(a); |
1136 | 1136 |
_next_out[tip] = a; |
1137 | 1137 |
tip = u; |
1138 | 1138 |
goto next_step; |
1139 | 1139 |
} |
0 comments (0 inline)