gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Fix and improve refine methods in CostScaling (#417)
0 1 0
default
1 file changed with 55 insertions and 34 deletions:
↑ Collapse diff ↑
Ignore white space 24 line context
... ...
@@ -1094,137 +1094,156 @@
1094 1094
        int k = std::min(_rank[u], r);
1095 1095
        if (k > 0) {
1096 1096
          _pi[u] -= _epsilon * k;
1097 1097
          _next_out[u] = _first_out[u];
1098 1098
        }
1099 1099
      }
1100 1100
    }
1101 1101

	
1102 1102
    /// Execute the algorithm performing augment and relabel operations
1103 1103
    void startAugment(int max_length) {
1104 1104
      // Paramters for heuristics
1105 1105
      const int EARLY_TERM_EPSILON_LIMIT = 1000;
1106
      const double GLOBAL_UPDATE_FACTOR = 3.0;
1107

	
1108
      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
1106
      const double GLOBAL_UPDATE_FACTOR = 1.0;
1107
      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
1109 1108
        (_res_node_num + _sup_node_num * _sup_node_num));
1110
      int next_update_limit = global_update_freq;
1111

	
1112
      int relabel_cnt = 0;
1109
      int next_global_update_limit = global_update_skip;
1113 1110

	
1114 1111
      // Perform cost scaling phases
1115
      std::vector<int> path;
1112
      IntVector path;
1113
      BoolVector path_arc(_res_arc_num, false);
1114
      int relabel_cnt = 0;
1116 1115
      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1117 1116
                                        1 : _epsilon / _alpha )
1118 1117
      {
1119 1118
        // Early termination heuristic
1120 1119
        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1121 1120
          if (earlyTermination()) break;
1122 1121
        }
1123 1122

	
1124 1123
        // Initialize current phase
1125 1124
        initPhase();
1126 1125

	
1127 1126
        // Perform partial augment and relabel operations
1128 1127
        while (true) {
1129 1128
          // Select an active node (FIFO selection)
1130 1129
          while (_active_nodes.size() > 0 &&
1131 1130
                 _excess[_active_nodes.front()] <= 0) {
1132 1131
            _active_nodes.pop_front();
1133 1132
          }
1134 1133
          if (_active_nodes.size() == 0) break;
1135 1134
          int start = _active_nodes.front();
1136 1135

	
1137 1136
          // Find an augmenting path from the start node
1138
          path.clear();
1139 1137
          int tip = start;
1140
          while (_excess[tip] >= 0 && int(path.size()) < max_length) {
1138
          while (int(path.size()) < max_length && _excess[tip] >= 0) {
1141 1139
            int u;
1142
            LargeCost min_red_cost, rc, pi_tip = _pi[tip];
1140
            LargeCost rc, min_red_cost = std::numeric_limits<LargeCost>::max();
1141
            LargeCost pi_tip = _pi[tip];
1143 1142
            int last_out = _first_out[tip+1];
1144 1143
            for (int a = _next_out[tip]; a != last_out; ++a) {
1145
              u = _target[a];
1146
              if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
1147
                path.push_back(a);
1148
                _next_out[tip] = a;
1149
                tip = u;
1150
                goto next_step;
1144
              if (_res_cap[a] > 0) {
1145
                u = _target[a];
1146
                rc = _cost[a] + pi_tip - _pi[u];
1147
                if (rc < 0) {
1148
                  path.push_back(a);
1149
                  _next_out[tip] = a;
1150
                  if (path_arc[a]) {
1151
                    goto augment;   // a cycle is found, stop path search
1152
                  }
1153
                  tip = u;
1154
                  path_arc[a] = true;
1155
                  goto next_step;
1156
                }
1157
                else if (rc < min_red_cost) {
1158
                  min_red_cost = rc;
1159
                }
1151 1160
              }
1152 1161
            }
1153 1162

	
1154 1163
            // Relabel tip node
1155
            min_red_cost = std::numeric_limits<LargeCost>::max();
1156 1164
            if (tip != start) {
1157 1165
              int ra = _reverse[path.back()];
1158
              min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
1166
              min_red_cost =
1167
                std::min(min_red_cost, _cost[ra] + pi_tip - _pi[_target[ra]]);
1159 1168
            }
1169
            last_out = _next_out[tip];
1160 1170
            for (int a = _first_out[tip]; a != last_out; ++a) {
1161
              rc = _cost[a] + pi_tip - _pi[_target[a]];
1162
              if (_res_cap[a] > 0 && rc < min_red_cost) {
1163
                min_red_cost = rc;
1171
              if (_res_cap[a] > 0) {
1172
                rc = _cost[a] + pi_tip - _pi[_target[a]];
1173
                if (rc < min_red_cost) {
1174
                  min_red_cost = rc;
1175
                }
1164 1176
              }
1165 1177
            }
1166 1178
            _pi[tip] -= min_red_cost + _epsilon;
1167 1179
            _next_out[tip] = _first_out[tip];
1168 1180
            ++relabel_cnt;
1169 1181

	
1170 1182
            // Step back
1171 1183
            if (tip != start) {
1172
              tip = _source[path.back()];
1184
              int pa = path.back();
1185
              path_arc[pa] = false;
1186
              tip = _source[pa];
1173 1187
              path.pop_back();
1174 1188
            }
1175 1189

	
1176 1190
          next_step: ;
1177 1191
          }
1178 1192

	
1179 1193
          // Augment along the found path (as much flow as possible)
1194
        augment:
1180 1195
          Value delta;
1181 1196
          int pa, u, v = start;
1182 1197
          for (int i = 0; i != int(path.size()); ++i) {
1183 1198
            pa = path[i];
1184 1199
            u = v;
1185 1200
            v = _target[pa];
1201
            path_arc[pa] = false;
1186 1202
            delta = std::min(_res_cap[pa], _excess[u]);
1187 1203
            _res_cap[pa] -= delta;
1188 1204
            _res_cap[_reverse[pa]] += delta;
1189 1205
            _excess[u] -= delta;
1190 1206
            _excess[v] += delta;
1191
            if (_excess[v] > 0 && _excess[v] <= delta)
1207
            if (_excess[v] > 0 && _excess[v] <= delta) {
1192 1208
              _active_nodes.push_back(v);
1209
            }
1193 1210
          }
1211
          path.clear();
1194 1212

	
1195 1213
          // Global update heuristic
1196
          if (relabel_cnt >= next_update_limit) {
1214
          if (relabel_cnt >= next_global_update_limit) {
1197 1215
            globalUpdate();
1198
            next_update_limit += global_update_freq;
1216
            next_global_update_limit += global_update_skip;
1199 1217
          }
1200 1218
        }
1219

	
1201 1220
      }
1221

	
1202 1222
    }
1203 1223

	
1204 1224
    /// Execute the algorithm performing push and relabel operations
1205 1225
    void startPush() {
1206 1226
      // Paramters for heuristics
1207 1227
      const int EARLY_TERM_EPSILON_LIMIT = 1000;
1208 1228
      const double GLOBAL_UPDATE_FACTOR = 2.0;
1209 1229

	
1210
      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
1230
      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
1211 1231
        (_res_node_num + _sup_node_num * _sup_node_num));
1212
      int next_update_limit = global_update_freq;
1213

	
1214
      int relabel_cnt = 0;
1232
      int next_global_update_limit = global_update_skip;
1215 1233

	
1216 1234
      // Perform cost scaling phases
1217 1235
      BoolVector hyper(_res_node_num, false);
1218 1236
      LargeCostVector hyper_cost(_res_node_num);
1237
      int relabel_cnt = 0;
1219 1238
      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1220 1239
                                        1 : _epsilon / _alpha )
1221 1240
      {
1222 1241
        // Early termination heuristic
1223 1242
        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1224 1243
          if (earlyTermination()) break;
1225 1244
        }
1226 1245

	
1227 1246
        // Initialize current phase
1228 1247
        initPhase();
1229 1248

	
1230 1249
        // Perform push and relabel operations
... ...
@@ -1284,49 +1303,51 @@
1284 1303
                  goto remove_nodes;
1285 1304
                }
1286 1305
              }
1287 1306
            }
1288 1307
            _next_out[n] = a;
1289 1308
          }
1290 1309

	
1291 1310
          // Relabel the node if it is still active (or hyper)
1292 1311
          if (_excess[n] > 0 || hyper[n]) {
1293 1312
             min_red_cost = hyper[n] ? -hyper_cost[n] :
1294 1313
               std::numeric_limits<LargeCost>::max();
1295 1314
            for (int a = _first_out[n]; a != last_out; ++a) {
1296
              rc = _cost[a] + pi_n - _pi[_target[a]];
1297
              if (_res_cap[a] > 0 && rc < min_red_cost) {
1298
                min_red_cost = rc;
1315
              if (_res_cap[a] > 0) {
1316
                rc = _cost[a] + pi_n - _pi[_target[a]];
1317
                if (rc < min_red_cost) {
1318
                  min_red_cost = rc;
1319
                }
1299 1320
              }
1300 1321
            }
1301 1322
            _pi[n] -= min_red_cost + _epsilon;
1302 1323
            _next_out[n] = _first_out[n];
1303 1324
            hyper[n] = false;
1304 1325
            ++relabel_cnt;
1305 1326
          }
1306 1327

	
1307 1328
          // Remove nodes that are not active nor hyper
1308 1329
        remove_nodes:
1309 1330
          while ( _active_nodes.size() > 0 &&
1310 1331
                  _excess[_active_nodes.front()] <= 0 &&
1311 1332
                  !hyper[_active_nodes.front()] ) {
1312 1333
            _active_nodes.pop_front();
1313 1334
          }
1314 1335

	
1315 1336
          // Global update heuristic
1316
          if (relabel_cnt >= next_update_limit) {
1337
          if (relabel_cnt >= next_global_update_limit) {
1317 1338
            globalUpdate();
1318 1339
            for (int u = 0; u != _res_node_num; ++u)
1319 1340
              hyper[u] = false;
1320
            next_update_limit += global_update_freq;
1341
            next_global_update_limit += global_update_skip;
1321 1342
          }
1322 1343
        }
1323 1344
      }
1324 1345
    }
1325 1346

	
1326 1347
  }; //class CostScaling
1327 1348

	
1328 1349
  ///@}
1329 1350

	
1330 1351
} //namespace lemon
1331 1352

	
1332 1353
#endif //LEMON_COST_SCALING_H
0 comments (0 inline)