573 } |
573 } |
574 _have_lower = false; |
574 _have_lower = false; |
575 return *this; |
575 return *this; |
576 } |
576 } |
577 |
577 |
578 /// \brief Reset all the parameters that have been given before. |
578 /// \brief Reset the internal data structures and all the parameters |
579 /// |
579 /// that have been given before. |
580 /// This function resets all the paramaters that have been given |
580 /// |
581 /// before using functions \ref lowerMap(), \ref upperMap(), |
581 /// This function resets the internal data structures and all the |
582 /// \ref costMap(), \ref supplyMap(), \ref stSupply(). |
582 /// paramaters that have been given before using functions \ref lowerMap(), |
583 /// |
583 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
584 /// It is useful for multiple run() calls. If this function is not |
584 /// |
585 /// used, all the parameters given before are kept for the next |
585 /// It is useful for multiple \ref run() calls. By default, all the given |
586 /// \ref run() call. |
586 /// parameters are kept for the next \ref run() call, unless |
587 /// However, the underlying digraph must not be modified after this |
587 /// \ref resetParams() or \ref reset() is used. |
588 /// class have been constructed, since it copies and extends the graph. |
588 /// If the underlying digraph was also modified after the construction |
|
589 /// of the class or the last \ref reset() call, then the \ref reset() |
|
590 /// function must be used, otherwise \ref resetParams() is sufficient. |
|
591 /// |
|
592 /// See \ref resetParams() for examples. |
|
593 /// |
589 /// \return <tt>(*this)</tt> |
594 /// \return <tt>(*this)</tt> |
|
595 /// |
|
596 /// \see resetParams(), run() |
590 CostScaling& reset() { |
597 CostScaling& reset() { |
591 // Resize vectors |
598 // Resize vectors |
592 _node_num = countNodes(_graph); |
599 _node_num = countNodes(_graph); |
593 _arc_num = countArcs(_graph); |
600 _arc_num = countArcs(_graph); |
594 _res_node_num = _node_num + 1; |
601 _res_node_num = _node_num + 1; |
888 _cost[a] = 0; |
895 _cost[a] = 0; |
889 _cost[ra] = 0; |
896 _cost[ra] = 0; |
890 } |
897 } |
891 } |
898 } |
892 |
899 |
893 return OPTIMAL; |
|
894 } |
|
895 |
|
896 // Execute the algorithm and transform the results |
|
897 void start(Method method) { |
|
898 // Maximum path length for partial augment |
|
899 const int MAX_PATH_LENGTH = 4; |
|
900 |
|
901 // Initialize data structures for buckets |
900 // Initialize data structures for buckets |
902 _max_rank = _alpha * _res_node_num; |
901 _max_rank = _alpha * _res_node_num; |
903 _buckets.resize(_max_rank); |
902 _buckets.resize(_max_rank); |
904 _bucket_next.resize(_res_node_num + 1); |
903 _bucket_next.resize(_res_node_num + 1); |
905 _bucket_prev.resize(_res_node_num + 1); |
904 _bucket_prev.resize(_res_node_num + 1); |
906 _rank.resize(_res_node_num + 1); |
905 _rank.resize(_res_node_num + 1); |
907 |
906 |
908 // Execute the algorithm |
907 return OPTIMAL; |
|
908 } |
|
909 |
|
910 // Execute the algorithm and transform the results |
|
911 void start(Method method) { |
|
912 const int MAX_PARTIAL_PATH_LENGTH = 4; |
|
913 |
909 switch (method) { |
914 switch (method) { |
910 case PUSH: |
915 case PUSH: |
911 startPush(); |
916 startPush(); |
912 break; |
917 break; |
913 case AUGMENT: |
918 case AUGMENT: |
914 startAugment(_res_node_num - 1); |
919 startAugment(_res_node_num - 1); |
915 break; |
920 break; |
916 case PARTIAL_AUGMENT: |
921 case PARTIAL_AUGMENT: |
917 startAugment(MAX_PATH_LENGTH); |
922 startAugment(MAX_PARTIAL_PATH_LENGTH); |
918 break; |
923 break; |
919 } |
924 } |
920 |
925 |
921 // Compute node potentials for the original costs |
926 // Compute node potentials for the original costs |
922 _arc_vec.clear(); |
927 _arc_vec.clear(); |
949 // Saturate arcs not satisfying the optimality condition |
954 // Saturate arcs not satisfying the optimality condition |
950 for (int u = 0; u != _res_node_num; ++u) { |
955 for (int u = 0; u != _res_node_num; ++u) { |
951 int last_out = _first_out[u+1]; |
956 int last_out = _first_out[u+1]; |
952 LargeCost pi_u = _pi[u]; |
957 LargeCost pi_u = _pi[u]; |
953 for (int a = _first_out[u]; a != last_out; ++a) { |
958 for (int a = _first_out[u]; a != last_out; ++a) { |
954 int v = _target[a]; |
959 Value delta = _res_cap[a]; |
955 if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) { |
960 if (delta > 0) { |
956 Value delta = _res_cap[a]; |
961 int v = _target[a]; |
957 _excess[u] -= delta; |
962 if (_cost[a] + pi_u - _pi[v] < 0) { |
958 _excess[v] += delta; |
963 _excess[u] -= delta; |
959 _res_cap[a] = 0; |
964 _excess[v] += delta; |
960 _res_cap[_reverse[a]] += delta; |
965 _res_cap[a] = 0; |
|
966 _res_cap[_reverse[a]] += delta; |
|
967 } |
961 } |
968 } |
962 } |
969 } |
963 } |
970 } |
964 |
971 |
965 // Find active nodes (i.e. nodes with positive excess) |
972 // Find active nodes (i.e. nodes with positive excess) |
999 return done; |
1006 return done; |
1000 } |
1007 } |
1001 |
1008 |
1002 // Global potential update heuristic |
1009 // Global potential update heuristic |
1003 void globalUpdate() { |
1010 void globalUpdate() { |
1004 int bucket_end = _root + 1; |
1011 const int bucket_end = _root + 1; |
1005 |
1012 |
1006 // Initialize buckets |
1013 // Initialize buckets |
1007 for (int r = 0; r != _max_rank; ++r) { |
1014 for (int r = 0; r != _max_rank; ++r) { |
1008 _buckets[r] = bucket_end; |
1015 _buckets[r] = bucket_end; |
1009 } |
1016 } |
1010 Value total_excess = 0; |
1017 Value total_excess = 0; |
|
1018 int b0 = bucket_end; |
1011 for (int i = 0; i != _res_node_num; ++i) { |
1019 for (int i = 0; i != _res_node_num; ++i) { |
1012 if (_excess[i] < 0) { |
1020 if (_excess[i] < 0) { |
1013 _rank[i] = 0; |
1021 _rank[i] = 0; |
1014 _bucket_next[i] = _buckets[0]; |
1022 _bucket_next[i] = b0; |
1015 _bucket_prev[_buckets[0]] = i; |
1023 _bucket_prev[b0] = i; |
1016 _buckets[0] = i; |
1024 b0 = i; |
1017 } else { |
1025 } else { |
1018 total_excess += _excess[i]; |
1026 total_excess += _excess[i]; |
1019 _rank[i] = _max_rank; |
1027 _rank[i] = _max_rank; |
1020 } |
1028 } |
1021 } |
1029 } |
1022 if (total_excess == 0) return; |
1030 if (total_excess == 0) return; |
|
1031 _buckets[0] = b0; |
1023 |
1032 |
1024 // Search the buckets |
1033 // Search the buckets |
1025 int r = 0; |
1034 int r = 0; |
1026 for ( ; r != _max_rank; ++r) { |
1035 for ( ; r != _max_rank; ++r) { |
1027 while (_buckets[r] != bucket_end) { |
1036 while (_buckets[r] != bucket_end) { |
1039 int old_rank_v = _rank[v]; |
1048 int old_rank_v = _rank[v]; |
1040 if (r < old_rank_v) { |
1049 if (r < old_rank_v) { |
1041 // Compute the new rank of v |
1050 // Compute the new rank of v |
1042 LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon; |
1051 LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon; |
1043 int new_rank_v = old_rank_v; |
1052 int new_rank_v = old_rank_v; |
1044 if (nrc < LargeCost(_max_rank)) |
1053 if (nrc < LargeCost(_max_rank)) { |
1045 new_rank_v = r + 1 + int(nrc); |
1054 new_rank_v = r + 1 + static_cast<int>(nrc); |
|
1055 } |
1046 |
1056 |
1047 // Change the rank of v |
1057 // Change the rank of v |
1048 if (new_rank_v < old_rank_v) { |
1058 if (new_rank_v < old_rank_v) { |
1049 _rank[v] = new_rank_v; |
1059 _rank[v] = new_rank_v; |
1050 _next_out[v] = _first_out[v]; |
1060 _next_out[v] = _first_out[v]; |
1052 // Remove v from its old bucket |
1062 // Remove v from its old bucket |
1053 if (old_rank_v < _max_rank) { |
1063 if (old_rank_v < _max_rank) { |
1054 if (_buckets[old_rank_v] == v) { |
1064 if (_buckets[old_rank_v] == v) { |
1055 _buckets[old_rank_v] = _bucket_next[v]; |
1065 _buckets[old_rank_v] = _bucket_next[v]; |
1056 } else { |
1066 } else { |
1057 _bucket_next[_bucket_prev[v]] = _bucket_next[v]; |
1067 int pv = _bucket_prev[v], nv = _bucket_next[v]; |
1058 _bucket_prev[_bucket_next[v]] = _bucket_prev[v]; |
1068 _bucket_next[pv] = nv; |
|
1069 _bucket_prev[nv] = pv; |
1059 } |
1070 } |
1060 } |
1071 } |
1061 |
1072 |
1062 // Insert v to its new bucket |
1073 // Insert v into its new bucket |
1063 _bucket_next[v] = _buckets[new_rank_v]; |
1074 int nv = _buckets[new_rank_v]; |
1064 _bucket_prev[_buckets[new_rank_v]] = v; |
1075 _bucket_next[v] = nv; |
|
1076 _bucket_prev[nv] = v; |
1065 _buckets[new_rank_v] = v; |
1077 _buckets[new_rank_v] = v; |
1066 } |
1078 } |
1067 } |
1079 } |
1068 } |
1080 } |
1069 } |
1081 } |