gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Improve the tree update process and a pivot rule (#391) and make some parts of the code clearer using better names
0 1 0
default
1 file changed with 139 insertions and 124 deletions:
↑ Collapse diff ↑
Ignore white space 12 line context
... ...
@@ -163,25 +163,28 @@
163 163

	
164 164
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
165 165

	
166 166
    typedef std::vector<int> IntVector;
167 167
    typedef std::vector<Value> ValueVector;
168 168
    typedef std::vector<Cost> CostVector;
169
    typedef std::vector<char> BoolVector;
170
    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
169
    typedef std::vector<signed char> CharVector;
170
    // Note: vector<signed char> is used instead of vector<ArcState> and 
171
    // vector<ArcDirection> for efficiency reasons
171 172

	
172 173
    // State constants for arcs
173 174
    enum ArcState {
174 175
      STATE_UPPER = -1,
175 176
      STATE_TREE  =  0,
176 177
      STATE_LOWER =  1
177 178
    };
178 179

	
179
    typedef std::vector<signed char> StateVector;
180
    // Note: vector<signed char> is used instead of vector<ArcState> for
181
    // efficiency reasons
180
    // Direction constants for tree arcs
181
    enum ArcDirection {
182
      DIR_DOWN = -1,
183
      DIR_UP   =  1
184
    };
182 185

	
183 186
  private:
184 187

	
185 188
    // Data related to the underlying digraph
186 189
    const GR &_graph;
187 190
    int _node_num;
... ...
@@ -214,21 +217,19 @@
214 217
    IntVector _parent;
215 218
    IntVector _pred;
216 219
    IntVector _thread;
217 220
    IntVector _rev_thread;
218 221
    IntVector _succ_num;
219 222
    IntVector _last_succ;
223
    CharVector _pred_dir;
224
    CharVector _state;
220 225
    IntVector _dirty_revs;
221
    BoolVector _forward;
222
    StateVector _state;
223 226
    int _root;
224 227

	
225 228
    // Temporary data used in the current pivot iteration
226 229
    int in_arc, join, u_in, v_in, u_out, v_out;
227
    int first, second, right, last;
228
    int stem, par_stem, new_stem;
229 230
    Value delta;
230 231

	
231 232
    const Value MAX;
232 233

	
233 234
  public:
234 235

	
... ...
@@ -247,13 +248,13 @@
247 248
    private:
248 249

	
249 250
      // References to the NetworkSimplex class
250 251
      const IntVector  &_source;
251 252
      const IntVector  &_target;
252 253
      const CostVector &_cost;
253
      const StateVector &_state;
254
      const CharVector &_state;
254 255
      const CostVector &_pi;
255 256
      int &_in_arc;
256 257
      int _search_arc_num;
257 258

	
258 259
      // Pivot rule data
259 260
      int _next_arc;
... ...
@@ -299,13 +300,13 @@
299 300
    private:
300 301

	
301 302
      // References to the NetworkSimplex class
302 303
      const IntVector  &_source;
303 304
      const IntVector  &_target;
304 305
      const CostVector &_cost;
305
      const StateVector &_state;
306
      const CharVector &_state;
306 307
      const CostVector &_pi;
307 308
      int &_in_arc;
308 309
      int _search_arc_num;
309 310

	
310 311
    public:
311 312

	
... ...
@@ -338,13 +339,13 @@
338 339
    private:
339 340

	
340 341
      // References to the NetworkSimplex class
341 342
      const IntVector  &_source;
342 343
      const IntVector  &_target;
343 344
      const CostVector &_cost;
344
      const StateVector &_state;
345
      const CharVector &_state;
345 346
      const CostVector &_pi;
346 347
      int &_in_arc;
347 348
      int _search_arc_num;
348 349

	
349 350
      // Pivot rule data
350 351
      int _block_size;
... ...
@@ -411,13 +412,13 @@
411 412
    private:
412 413

	
413 414
      // References to the NetworkSimplex class
414 415
      const IntVector  &_source;
415 416
      const IntVector  &_target;
416 417
      const CostVector &_cost;
417
      const StateVector &_state;
418
      const CharVector &_state;
418 419
      const CostVector &_pi;
419 420
      int &_in_arc;
420 421
      int _search_arc_num;
421 422

	
422 423
      // Pivot rule data
423 424
      IntVector _candidates;
... ...
@@ -514,13 +515,13 @@
514 515
    private:
515 516

	
516 517
      // References to the NetworkSimplex class
517 518
      const IntVector  &_source;
518 519
      const IntVector  &_target;
519 520
      const CostVector &_cost;
520
      const StateVector &_state;
521
      const CharVector &_state;
521 522
      const CostVector &_pi;
522 523
      int &_in_arc;
523 524
      int _search_arc_num;
524 525

	
525 526
      // Pivot rule data
526 527
      int _block_size, _head_length, _curr_length;
... ...
@@ -567,29 +568,31 @@
567 568
      }
568 569

	
569 570
      // Find next entering arc
570 571
      bool findEnteringArc() {
571 572
        // Check the current candidate list
572 573
        int e;
574
        Cost c;
573 575
        for (int i = 0; i != _curr_length; ++i) {
574 576
          e = _candidates[i];
575
          _cand_cost[e] = _state[e] *
576
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
577
          if (_cand_cost[e] >= 0) {
577
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
578
          if (c < 0) {
579
            _cand_cost[e] = c;
580
          } else {
578 581
            _candidates[i--] = _candidates[--_curr_length];
579 582
          }
580 583
        }
581 584

	
582 585
        // Extend the list
583 586
        int cnt = _block_size;
584 587
        int limit = _head_length;
585 588

	
586 589
        for (e = _next_arc; e != _search_arc_num; ++e) {
587
          _cand_cost[e] = _state[e] *
588
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
589
          if (_cand_cost[e] < 0) {
590
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
591
          if (c < 0) {
592
            _cand_cost[e] = c;
590 593
            _candidates[_curr_length++] = e;
591 594
          }
592 595
          if (--cnt == 0) {
593 596
            if (_curr_length > limit) goto search_end;
594 597
            limit = 0;
595 598
            cnt = _block_size;
... ...
@@ -910,13 +913,13 @@
910 913
      _supply.resize(all_node_num);
911 914
      _flow.resize(max_arc_num);
912 915
      _pi.resize(all_node_num);
913 916

	
914 917
      _parent.resize(all_node_num);
915 918
      _pred.resize(all_node_num);
916
      _forward.resize(all_node_num);
919
      _pred_dir.resize(all_node_num);
917 920
      _thread.resize(all_node_num);
918 921
      _rev_thread.resize(all_node_num);
919 922
      _succ_num.resize(all_node_num);
920 923
      _last_succ.resize(all_node_num);
921 924
      _state.resize(max_arc_num);
922 925

	
... ...
@@ -1113,20 +1116,20 @@
1113 1116
          _rev_thread[u + 1] = u;
1114 1117
          _succ_num[u] = 1;
1115 1118
          _last_succ[u] = u;
1116 1119
          _cap[e] = INF;
1117 1120
          _state[e] = STATE_TREE;
1118 1121
          if (_supply[u] >= 0) {
1119
            _forward[u] = true;
1122
            _pred_dir[u] = DIR_UP;
1120 1123
            _pi[u] = 0;
1121 1124
            _source[e] = u;
1122 1125
            _target[e] = _root;
1123 1126
            _flow[e] = _supply[u];
1124 1127
            _cost[e] = 0;
1125 1128
          } else {
1126
            _forward[u] = false;
1129
            _pred_dir[u] = DIR_DOWN;
1127 1130
            _pi[u] = ART_COST;
1128 1131
            _source[e] = _root;
1129 1132
            _target[e] = u;
1130 1133
            _flow[e] = -_supply[u];
1131 1134
            _cost[e] = ART_COST;
1132 1135
          }
... ...
@@ -1140,23 +1143,23 @@
1140 1143
          _parent[u] = _root;
1141 1144
          _thread[u] = u + 1;
1142 1145
          _rev_thread[u + 1] = u;
1143 1146
          _succ_num[u] = 1;
1144 1147
          _last_succ[u] = u;
1145 1148
          if (_supply[u] >= 0) {
1146
            _forward[u] = true;
1149
            _pred_dir[u] = DIR_UP;
1147 1150
            _pi[u] = 0;
1148 1151
            _pred[u] = e;
1149 1152
            _source[e] = u;
1150 1153
            _target[e] = _root;
1151 1154
            _cap[e] = INF;
1152 1155
            _flow[e] = _supply[u];
1153 1156
            _cost[e] = 0;
1154 1157
            _state[e] = STATE_TREE;
1155 1158
          } else {
1156
            _forward[u] = false;
1159
            _pred_dir[u] = DIR_DOWN;
1157 1160
            _pi[u] = ART_COST;
1158 1161
            _pred[u] = f;
1159 1162
            _source[f] = _root;
1160 1163
            _target[f] = u;
1161 1164
            _cap[f] = INF;
1162 1165
            _flow[f] = -_supply[u];
... ...
@@ -1181,23 +1184,23 @@
1181 1184
          _parent[u] = _root;
1182 1185
          _thread[u] = u + 1;
1183 1186
          _rev_thread[u + 1] = u;
1184 1187
          _succ_num[u] = 1;
1185 1188
          _last_succ[u] = u;
1186 1189
          if (_supply[u] <= 0) {
1187
            _forward[u] = false;
1190
            _pred_dir[u] = DIR_DOWN;
1188 1191
            _pi[u] = 0;
1189 1192
            _pred[u] = e;
1190 1193
            _source[e] = _root;
1191 1194
            _target[e] = u;
1192 1195
            _cap[e] = INF;
1193 1196
            _flow[e] = -_supply[u];
1194 1197
            _cost[e] = 0;
1195 1198
            _state[e] = STATE_TREE;
1196 1199
          } else {
1197
            _forward[u] = true;
1200
            _pred_dir[u] = DIR_UP;
1198 1201
            _pi[u] = -ART_COST;
1199 1202
            _pred[u] = f;
1200 1203
            _source[f] = u;
1201 1204
            _target[f] = _root;
1202 1205
            _cap[f] = INF;
1203 1206
            _flow[f] = _supply[u];
... ...
@@ -1234,40 +1237,48 @@
1234 1237

	
1235 1238
    // Find the leaving arc of the cycle and returns true if the
1236 1239
    // leaving arc is not the same as the entering arc
1237 1240
    bool findLeavingArc() {
1238 1241
      // Initialize first and second nodes according to the direction
1239 1242
      // of the cycle
1243
      int first, second;
1240 1244
      if (_state[in_arc] == STATE_LOWER) {
1241 1245
        first  = _source[in_arc];
1242 1246
        second = _target[in_arc];
1243 1247
      } else {
1244 1248
        first  = _target[in_arc];
1245 1249
        second = _source[in_arc];
1246 1250
      }
1247 1251
      delta = _cap[in_arc];
1248 1252
      int result = 0;
1249
      Value d;
1253
      Value c, d;
1250 1254
      int e;
1251 1255

	
1252
      // Search the cycle along the path form the first node to the root
1256
      // Search the cycle form the first node to the join node
1253 1257
      for (int u = first; u != join; u = _parent[u]) {
1254 1258
        e = _pred[u];
1255
        d = _forward[u] ?
1256
          _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
1259
        d = _flow[e];
1260
        if (_pred_dir[u] == DIR_DOWN) {
1261
          c = _cap[e];
1262
          d = c >= MAX ? INF : c - d;
1263
        }
1257 1264
        if (d < delta) {
1258 1265
          delta = d;
1259 1266
          u_out = u;
1260 1267
          result = 1;
1261 1268
        }
1262 1269
      }
1263
      // Search the cycle along the path form the second node to the root
1270

	
1271
      // Search the cycle form the second node to the join node
1264 1272
      for (int u = second; u != join; u = _parent[u]) {
1265 1273
        e = _pred[u];
1266
        d = _forward[u] ?
1267
          (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
1274
        d = _flow[e];
1275
        if (_pred_dir[u] == DIR_UP) {
1276
          c = _cap[e];
1277
          d = c >= MAX ? INF : c - d;
1278
        }
1268 1279
        if (d <= delta) {
1269 1280
          delta = d;
1270 1281
          u_out = u;
1271 1282
          result = 2;
1272 1283
        }
1273 1284
      }
... ...
@@ -1286,16 +1297,16 @@
1286 1297
    void changeFlow(bool change) {
1287 1298
      // Augment along the cycle
1288 1299
      if (delta > 0) {
1289 1300
        Value val = _state[in_arc] * delta;
1290 1301
        _flow[in_arc] += val;
1291 1302
        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
1292
          _flow[_pred[u]] += _forward[u] ? -val : val;
1303
          _flow[_pred[u]] -= _pred_dir[u] * val;
1293 1304
        }
1294 1305
        for (int u = _target[in_arc]; u != join; u = _parent[u]) {
1295
          _flow[_pred[u]] += _forward[u] ? val : -val;
1306
          _flow[_pred[u]] += _pred_dir[u] * val;
1296 1307
        }
1297 1308
      }
1298 1309
      // Update the state of the entering and leaving arcs
1299 1310
      if (change) {
1300 1311
        _state[in_arc] = STATE_TREE;
1301 1312
        _state[_pred[u_out]] =
... ...
@@ -1304,136 +1315,140 @@
1304 1315
        _state[in_arc] = -_state[in_arc];
1305 1316
      }
1306 1317
    }
1307 1318

	
1308 1319
    // Update the tree structure
1309 1320
    void updateTreeStructure() {
1310
      int u, w;
1311 1321
      int old_rev_thread = _rev_thread[u_out];
1312 1322
      int old_succ_num = _succ_num[u_out];
1313 1323
      int old_last_succ = _last_succ[u_out];
1314 1324
      v_out = _parent[u_out];
1315 1325

	
1316
      u = _last_succ[u_in];  // the last successor of u_in
1317
      right = _thread[u];    // the node after it
1326
      // Check if u_in and u_out coincide
1327
      if (u_in == u_out) {
1328
        // Update _parent, _pred, _pred_dir
1329
        _parent[u_in] = v_in;
1330
        _pred[u_in] = in_arc;
1331
        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
1318 1332

	
1319
      // Handle the case when old_rev_thread equals to v_in
1320
      // (it also means that join and v_out coincide)
1321
      if (old_rev_thread == v_in) {
1322
        last = _thread[_last_succ[u_out]];
1333
        // Update _thread and _rev_thread
1334
        if (_thread[v_in] != u_out) {
1335
          int after = _thread[old_last_succ];
1336
          _thread[old_rev_thread] = after;
1337
          _rev_thread[after] = old_rev_thread;
1338
          after = _thread[v_in];
1339
          _thread[v_in] = u_out;
1340
          _rev_thread[u_out] = v_in;
1341
          _thread[old_last_succ] = after;
1342
          _rev_thread[after] = old_last_succ;
1343
        }
1323 1344
      } else {
1324
        last = _thread[v_in];
1325
      }
1345
        // Handle the case when old_rev_thread equals to v_in
1346
        // (it also means that join and v_out coincide)
1347
        int thread_continue = old_rev_thread == v_in ?
1348
          _thread[old_last_succ] : _thread[v_in];
1326 1349

	
1327
      // Update _thread and _parent along the stem nodes (i.e. the nodes
1328
      // between u_in and u_out, whose parent have to be changed)
1329
      _thread[v_in] = stem = u_in;
1330
      _dirty_revs.clear();
1331
      _dirty_revs.push_back(v_in);
1332
      par_stem = v_in;
1333
      while (stem != u_out) {
1334
        // Insert the next stem node into the thread list
1335
        new_stem = _parent[stem];
1336
        _thread[u] = new_stem;
1337
        _dirty_revs.push_back(u);
1350
        // Update _thread and _parent along the stem nodes (i.e. the nodes
1351
        // between u_in and u_out, whose parent have to be changed)
1352
        int stem = u_in;              // the current stem node
1353
        int par_stem = v_in;          // the new parent of stem
1354
        int next_stem;                // the next stem node
1355
        int last = _last_succ[u_in];  // the last successor of stem
1356
        int before, after = _thread[last];
1357
        _thread[v_in] = u_in;
1358
        _dirty_revs.clear();
1359
        _dirty_revs.push_back(v_in);
1360
        while (stem != u_out) {
1361
          // Insert the next stem node into the thread list
1362
          next_stem = _parent[stem];
1363
          _thread[last] = next_stem;
1364
          _dirty_revs.push_back(last);
1338 1365

	
1339
        // Remove the subtree of stem from the thread list
1340
        w = _rev_thread[stem];
1341
        _thread[w] = right;
1342
        _rev_thread[right] = w;
1366
          // Remove the subtree of stem from the thread list
1367
          before = _rev_thread[stem];
1368
          _thread[before] = after;
1369
          _rev_thread[after] = before;
1343 1370

	
1344
        // Change the parent node and shift stem nodes
1345
        _parent[stem] = par_stem;
1346
        par_stem = stem;
1347
        stem = new_stem;
1371
          // Change the parent node and shift stem nodes
1372
          _parent[stem] = par_stem;
1373
          par_stem = stem;
1374
          stem = next_stem;
1348 1375

	
1349
        // Update u and right
1350
        u = _last_succ[stem] == _last_succ[par_stem] ?
1351
          _rev_thread[par_stem] : _last_succ[stem];
1352
        right = _thread[u];
1353
      }
1354
      _parent[u_out] = par_stem;
1355
      _thread[u] = last;
1356
      _rev_thread[last] = u;
1357
      _last_succ[u_out] = u;
1376
          // Update last and after
1377
          last = _last_succ[stem] == _last_succ[par_stem] ?
1378
            _rev_thread[par_stem] : _last_succ[stem];
1379
          after = _thread[last];
1380
        }
1381
        _parent[u_out] = par_stem;
1382
        _thread[last] = thread_continue;
1383
        _rev_thread[thread_continue] = last;
1384
        _last_succ[u_out] = last;
1358 1385

	
1359
      // Remove the subtree of u_out from the thread list except for
1360
      // the case when old_rev_thread equals to v_in
1361
      // (it also means that join and v_out coincide)
1362
      if (old_rev_thread != v_in) {
1363
        _thread[old_rev_thread] = right;
1364
        _rev_thread[right] = old_rev_thread;
1365
      }
1386
        // Remove the subtree of u_out from the thread list except for
1387
        // the case when old_rev_thread equals to v_in
1388
        if (old_rev_thread != v_in) {
1389
          _thread[old_rev_thread] = after;
1390
          _rev_thread[after] = old_rev_thread;
1391
        }
1366 1392

	
1367
      // Update _rev_thread using the new _thread values
1368
      for (int i = 0; i != int(_dirty_revs.size()); ++i) {
1369
        u = _dirty_revs[i];
1370
        _rev_thread[_thread[u]] = u;
1371
      }
1393
        // Update _rev_thread using the new _thread values
1394
        for (int i = 0; i != int(_dirty_revs.size()); ++i) {
1395
          int u = _dirty_revs[i];
1396
          _rev_thread[_thread[u]] = u;
1397
        }
1372 1398

	
1373
      // Update _pred, _forward, _last_succ and _succ_num for the
1374
      // stem nodes from u_out to u_in
1375
      int tmp_sc = 0, tmp_ls = _last_succ[u_out];
1376
      u = u_out;
1377
      while (u != u_in) {
1378
        w = _parent[u];
1379
        _pred[u] = _pred[w];
1380
        _forward[u] = !_forward[w];
1381
        tmp_sc += _succ_num[u] - _succ_num[w];
1382
        _succ_num[u] = tmp_sc;
1383
        _last_succ[w] = tmp_ls;
1384
        u = w;
1385
      }
1386
      _pred[u_in] = in_arc;
1387
      _forward[u_in] = (u_in == _source[in_arc]);
1388
      _succ_num[u_in] = old_succ_num;
1389

	
1390
      // Set limits for updating _last_succ form v_in and v_out
1391
      // towards the root
1392
      int up_limit_in = -1;
1393
      int up_limit_out = -1;
1394
      if (_last_succ[join] == v_in) {
1395
        up_limit_out = join;
1396
      } else {
1397
        up_limit_in = join;
1399
        // Update _pred, _pred_dir, _last_succ and _succ_num for the
1400
        // stem nodes from u_out to u_in
1401
        int tmp_sc = 0, tmp_ls = _last_succ[u_out];
1402
        for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) {
1403
          _pred[u] = _pred[p];
1404
          _pred_dir[u] = -_pred_dir[p];
1405
          tmp_sc += _succ_num[u] - _succ_num[p];
1406
          _succ_num[u] = tmp_sc;
1407
          _last_succ[p] = tmp_ls;
1408
        }
1409
        _pred[u_in] = in_arc;
1410
        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
1411
        _succ_num[u_in] = old_succ_num;
1398 1412
      }
1399 1413

	
1400 1414
      // Update _last_succ from v_in towards the root
1401
      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
1402
           u = _parent[u]) {
1403
        _last_succ[u] = _last_succ[u_out];
1415
      int up_limit_out = _last_succ[join] == v_in ? join : -1;
1416
      int last_succ_out = _last_succ[u_out];
1417
      for (int u = v_in; u != -1 && _last_succ[u] == v_in; u = _parent[u]) {
1418
        _last_succ[u] = last_succ_out;
1404 1419
      }
1420

	
1405 1421
      // Update _last_succ from v_out towards the root
1406 1422
      if (join != old_rev_thread && v_in != old_rev_thread) {
1407
        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1423
        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1408 1424
             u = _parent[u]) {
1409 1425
          _last_succ[u] = old_rev_thread;
1410 1426
        }
1411
      } else {
1412
        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1427
      }
1428
      else if (last_succ_out != old_last_succ) {
1429
        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1413 1430
             u = _parent[u]) {
1414
          _last_succ[u] = _last_succ[u_out];
1431
          _last_succ[u] = last_succ_out;
1415 1432
        }
1416 1433
      }
1417 1434

	
1418 1435
      // Update _succ_num from v_in to join
1419
      for (u = v_in; u != join; u = _parent[u]) {
1436
      for (int u = v_in; u != join; u = _parent[u]) {
1420 1437
        _succ_num[u] += old_succ_num;
1421 1438
      }
1422 1439
      // Update _succ_num from v_out to join
1423
      for (u = v_out; u != join; u = _parent[u]) {
1440
      for (int u = v_out; u != join; u = _parent[u]) {
1424 1441
        _succ_num[u] -= old_succ_num;
1425 1442
      }
1426 1443
    }
1427 1444

	
1428
    // Update potentials
1445
    // Update potentials in the subtree that has been moved
1429 1446
    void updatePotential() {
1430
      Cost sigma = _forward[u_in] ?
1431
        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1432
        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1433
      // Update potentials in the subtree, which has been moved
1447
      Cost sigma = _pi[v_in] - _pi[u_in] -
1448
                   _pred_dir[u_in] * _cost[in_arc];
1434 1449
      int end = _thread[_last_succ[u_in]];
1435 1450
      for (int u = u_in; u != end; u = _thread[u]) {
1436 1451
        _pi[u] += sigma;
1437 1452
      }
1438 1453
    }
1439 1454

	
0 comments (0 inline)