gravatar
alpar (Alpar Juttner)
alpar@cs.elte.hu
Merge #417
0 1 0
merge default
1 file changed with 2 insertions and 2 deletions:
↑ Collapse diff ↑
Ignore white space 96 line context
... ...
@@ -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 = std::numeric_limits<int>::max()) {
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)