gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Exploit that the standard maps are reference maps (#190)
0 9 0
default
9 files changed with 345 insertions and 345 deletions:
↑ Collapse diff ↑
Show white space 96 line context
... ...
@@ -408,223 +408,223 @@
408 408
        _local_level = false;
409 409
      }
410 410
      _level = &elevator;
411 411
      return *this;
412 412
    }
413 413

	
414 414
    /// \brief Returns a const reference to the elevator.
415 415
    ///
416 416
    /// Returns a const reference to the elevator.
417 417
    ///
418 418
    /// \pre Either \ref run() or \ref init() must be called before
419 419
    /// using this function.
420 420
    const Elevator& elevator() const {
421 421
      return *_level;
422 422
    }
423 423

	
424 424
    /// \brief Sets the tolerance used by algorithm.
425 425
    ///
426 426
    /// Sets the tolerance used by algorithm.
427 427
    Circulation& tolerance(const Tolerance& tolerance) const {
428 428
      _tol = tolerance;
429 429
      return *this;
430 430
    }
431 431

	
432 432
    /// \brief Returns a const reference to the tolerance.
433 433
    ///
434 434
    /// Returns a const reference to the tolerance.
435 435
    const Tolerance& tolerance() const {
436 436
      return tolerance;
437 437
    }
438 438

	
439 439
    /// \name Execution Control
440 440
    /// The simplest way to execute the algorithm is to call \ref run().\n
441 441
    /// If you need more control on the initial solution or the execution,
442 442
    /// first you have to call one of the \ref init() functions, then
443 443
    /// the \ref start() function.
444 444

	
445 445
    ///@{
446 446

	
447 447
    /// Initializes the internal data structures.
448 448

	
449 449
    /// Initializes the internal data structures and sets all flow values
450 450
    /// to the lower bound.
451 451
    void init()
452 452
    {
453 453
      createStructures();
454 454

	
455 455
      for(NodeIt n(_g);n!=INVALID;++n) {
456
        _excess->set(n, (*_delta)[n]);
456
        (*_excess)[n] = (*_delta)[n];
457 457
      }
458 458

	
459 459
      for (ArcIt e(_g);e!=INVALID;++e) {
460 460
        _flow->set(e, (*_lo)[e]);
461
        _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
462
        _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
461
        (*_excess)[_g.target(e)] += (*_flow)[e];
462
        (*_excess)[_g.source(e)] -= (*_flow)[e];
463 463
      }
464 464

	
465 465
      // global relabeling tested, but in general case it provides
466 466
      // worse performance for random digraphs
467 467
      _level->initStart();
468 468
      for(NodeIt n(_g);n!=INVALID;++n)
469 469
        _level->initAddItem(n);
470 470
      _level->initFinish();
471 471
      for(NodeIt n(_g);n!=INVALID;++n)
472 472
        if(_tol.positive((*_excess)[n]))
473 473
          _level->activate(n);
474 474
    }
475 475

	
476 476
    /// Initializes the internal data structures using a greedy approach.
477 477

	
478 478
    /// Initializes the internal data structures using a greedy approach
479 479
    /// to construct the initial solution.
480 480
    void greedyInit()
481 481
    {
482 482
      createStructures();
483 483

	
484 484
      for(NodeIt n(_g);n!=INVALID;++n) {
485
        _excess->set(n, (*_delta)[n]);
485
        (*_excess)[n] = (*_delta)[n];
486 486
      }
487 487

	
488 488
      for (ArcIt e(_g);e!=INVALID;++e) {
489 489
        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
490 490
          _flow->set(e, (*_up)[e]);
491
          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
492
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
491
          (*_excess)[_g.target(e)] += (*_up)[e];
492
          (*_excess)[_g.source(e)] -= (*_up)[e];
493 493
        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
494 494
          _flow->set(e, (*_lo)[e]);
495
          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
496
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
495
          (*_excess)[_g.target(e)] += (*_lo)[e];
496
          (*_excess)[_g.source(e)] -= (*_lo)[e];
497 497
        } else {
498 498
          Value fc = -(*_excess)[_g.target(e)];
499 499
          _flow->set(e, fc);
500
          _excess->set(_g.target(e), 0);
501
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
500
          (*_excess)[_g.target(e)] = 0;
501
          (*_excess)[_g.source(e)] -= fc;
502 502
        }
503 503
      }
504 504

	
505 505
      _level->initStart();
506 506
      for(NodeIt n(_g);n!=INVALID;++n)
507 507
        _level->initAddItem(n);
508 508
      _level->initFinish();
509 509
      for(NodeIt n(_g);n!=INVALID;++n)
510 510
        if(_tol.positive((*_excess)[n]))
511 511
          _level->activate(n);
512 512
    }
513 513

	
514 514
    ///Executes the algorithm
515 515

	
516 516
    ///This function executes the algorithm.
517 517
    ///
518 518
    ///\return \c true if a feasible circulation is found.
519 519
    ///
520 520
    ///\sa barrier()
521 521
    ///\sa barrierMap()
522 522
    bool start()
523 523
    {
524 524

	
525 525
      Node act;
526 526
      Node bact=INVALID;
527 527
      Node last_activated=INVALID;
528 528
      while((act=_level->highestActive())!=INVALID) {
529 529
        int actlevel=(*_level)[act];
530 530
        int mlevel=_node_num;
531 531
        Value exc=(*_excess)[act];
532 532

	
533 533
        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
534 534
          Node v = _g.target(e);
535 535
          Value fc=(*_up)[e]-(*_flow)[e];
536 536
          if(!_tol.positive(fc)) continue;
537 537
          if((*_level)[v]<actlevel) {
538 538
            if(!_tol.less(fc, exc)) {
539 539
              _flow->set(e, (*_flow)[e] + exc);
540
              _excess->set(v, (*_excess)[v] + exc);
540
              (*_excess)[v] += exc;
541 541
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
542 542
                _level->activate(v);
543
              _excess->set(act,0);
543
              (*_excess)[act] = 0;
544 544
              _level->deactivate(act);
545 545
              goto next_l;
546 546
            }
547 547
            else {
548 548
              _flow->set(e, (*_up)[e]);
549
              _excess->set(v, (*_excess)[v] + fc);
549
              (*_excess)[v] += fc;
550 550
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
551 551
                _level->activate(v);
552 552
              exc-=fc;
553 553
            }
554 554
          }
555 555
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
556 556
        }
557 557
        for(InArcIt e(_g,act);e!=INVALID; ++e) {
558 558
          Node v = _g.source(e);
559 559
          Value fc=(*_flow)[e]-(*_lo)[e];
560 560
          if(!_tol.positive(fc)) continue;
561 561
          if((*_level)[v]<actlevel) {
562 562
            if(!_tol.less(fc, exc)) {
563 563
              _flow->set(e, (*_flow)[e] - exc);
564
              _excess->set(v, (*_excess)[v] + exc);
564
              (*_excess)[v] += exc;
565 565
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
566 566
                _level->activate(v);
567
              _excess->set(act,0);
567
              (*_excess)[act] = 0;
568 568
              _level->deactivate(act);
569 569
              goto next_l;
570 570
            }
571 571
            else {
572 572
              _flow->set(e, (*_lo)[e]);
573
              _excess->set(v, (*_excess)[v] + fc);
573
              (*_excess)[v] += fc;
574 574
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
575 575
                _level->activate(v);
576 576
              exc-=fc;
577 577
            }
578 578
          }
579 579
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
580 580
        }
581 581

	
582
        _excess->set(act, exc);
582
        (*_excess)[act] = exc;
583 583
        if(!_tol.positive(exc)) _level->deactivate(act);
584 584
        else if(mlevel==_node_num) {
585 585
          _level->liftHighestActiveToTop();
586 586
          _el = _node_num;
587 587
          return false;
588 588
        }
589 589
        else {
590 590
          _level->liftHighestActive(mlevel+1);
591 591
          if(_level->onLevel(actlevel)==0) {
592 592
            _el = actlevel;
593 593
            return false;
594 594
          }
595 595
        }
596 596
      next_l:
597 597
        ;
598 598
      }
599 599
      return true;
600 600
    }
601 601

	
602 602
    /// Runs the algorithm.
603 603

	
604 604
    /// This function runs the algorithm.
605 605
    ///
606 606
    /// \return \c true if a feasible circulation is found.
607 607
    ///
608 608
    /// \note Apart from the return value, c.run() is just a shortcut of
609 609
    /// the following code.
610 610
    /// \code
611 611
    ///   c.greedyInit();
612 612
    ///   c.start();
613 613
    /// \endcode
614 614
    bool run() {
615 615
      greedyInit();
616 616
      return start();
617 617
    }
618 618

	
619 619
    /// @}
620 620

	
621 621
    /// \name Query Functions
622 622
    /// The results of the circulation algorithm can be obtained using
623 623
    /// these functions.\n
624 624
    /// Either \ref run() or \ref start() should be called before
625 625
    /// using them.
626 626

	
627 627
    ///@{
628 628

	
629 629
    /// \brief Returns the flow on the given arc.
630 630
    ///
Show white space 96 line context
... ...
@@ -1270,273 +1270,273 @@
1270 1270
      bool operator()(Arc a,Arc b) const
1271 1271
      {
1272 1272
        return g.target(a)<g.target(b);
1273 1273
      }
1274 1274
    };
1275 1275

	
1276 1276
  public:
1277 1277

	
1278 1278
    ///Constructor
1279 1279

	
1280 1280
    ///Constructor.
1281 1281
    ///
1282 1282
    ///It builds up the search database.
1283 1283
    DynArcLookUp(const Digraph &g)
1284 1284
      : _g(g),_head(g),_parent(g),_left(g),_right(g)
1285 1285
    {
1286 1286
      Parent::attach(_g.notifier(typename Digraph::Arc()));
1287 1287
      refresh();
1288 1288
    }
1289 1289

	
1290 1290
  protected:
1291 1291

	
1292 1292
    virtual void add(const Arc& arc) {
1293 1293
      insert(arc);
1294 1294
    }
1295 1295

	
1296 1296
    virtual void add(const std::vector<Arc>& arcs) {
1297 1297
      for (int i = 0; i < int(arcs.size()); ++i) {
1298 1298
        insert(arcs[i]);
1299 1299
      }
1300 1300
    }
1301 1301

	
1302 1302
    virtual void erase(const Arc& arc) {
1303 1303
      remove(arc);
1304 1304
    }
1305 1305

	
1306 1306
    virtual void erase(const std::vector<Arc>& arcs) {
1307 1307
      for (int i = 0; i < int(arcs.size()); ++i) {
1308 1308
        remove(arcs[i]);
1309 1309
      }
1310 1310
    }
1311 1311

	
1312 1312
    virtual void build() {
1313 1313
      refresh();
1314 1314
    }
1315 1315

	
1316 1316
    virtual void clear() {
1317 1317
      for(NodeIt n(_g);n!=INVALID;++n) {
1318
        _head.set(n, INVALID);
1318
        _head[n] = INVALID;
1319 1319
      }
1320 1320
    }
1321 1321

	
1322 1322
    void insert(Arc arc) {
1323 1323
      Node s = _g.source(arc);
1324 1324
      Node t = _g.target(arc);
1325
      _left.set(arc, INVALID);
1326
      _right.set(arc, INVALID);
1325
      _left[arc] = INVALID;
1326
      _right[arc] = INVALID;
1327 1327

	
1328 1328
      Arc e = _head[s];
1329 1329
      if (e == INVALID) {
1330
        _head.set(s, arc);
1331
        _parent.set(arc, INVALID);
1330
        _head[s] = arc;
1331
        _parent[arc] = INVALID;
1332 1332
        return;
1333 1333
      }
1334 1334
      while (true) {
1335 1335
        if (t < _g.target(e)) {
1336 1336
          if (_left[e] == INVALID) {
1337
            _left.set(e, arc);
1338
            _parent.set(arc, e);
1337
            _left[e] = arc;
1338
            _parent[arc] = e;
1339 1339
            splay(arc);
1340 1340
            return;
1341 1341
          } else {
1342 1342
            e = _left[e];
1343 1343
          }
1344 1344
        } else {
1345 1345
          if (_right[e] == INVALID) {
1346
            _right.set(e, arc);
1347
            _parent.set(arc, e);
1346
            _right[e] = arc;
1347
            _parent[arc] = e;
1348 1348
            splay(arc);
1349 1349
            return;
1350 1350
          } else {
1351 1351
            e = _right[e];
1352 1352
          }
1353 1353
        }
1354 1354
      }
1355 1355
    }
1356 1356

	
1357 1357
    void remove(Arc arc) {
1358 1358
      if (_left[arc] == INVALID) {
1359 1359
        if (_right[arc] != INVALID) {
1360
          _parent.set(_right[arc], _parent[arc]);
1360
          _parent[_right[arc]] = _parent[arc];
1361 1361
        }
1362 1362
        if (_parent[arc] != INVALID) {
1363 1363
          if (_left[_parent[arc]] == arc) {
1364
            _left.set(_parent[arc], _right[arc]);
1364
            _left[_parent[arc]] = _right[arc];
1365 1365
          } else {
1366
            _right.set(_parent[arc], _right[arc]);
1366
            _right[_parent[arc]] = _right[arc];
1367 1367
          }
1368 1368
        } else {
1369
          _head.set(_g.source(arc), _right[arc]);
1369
          _head[_g.source(arc)] = _right[arc];
1370 1370
        }
1371 1371
      } else if (_right[arc] == INVALID) {
1372
        _parent.set(_left[arc], _parent[arc]);
1372
        _parent[_left[arc]] = _parent[arc];
1373 1373
        if (_parent[arc] != INVALID) {
1374 1374
          if (_left[_parent[arc]] == arc) {
1375
            _left.set(_parent[arc], _left[arc]);
1375
            _left[_parent[arc]] = _left[arc];
1376 1376
          } else {
1377
            _right.set(_parent[arc], _left[arc]);
1377
            _right[_parent[arc]] = _left[arc];
1378 1378
          }
1379 1379
        } else {
1380
          _head.set(_g.source(arc), _left[arc]);
1380
          _head[_g.source(arc)] = _left[arc];
1381 1381
        }
1382 1382
      } else {
1383 1383
        Arc e = _left[arc];
1384 1384
        if (_right[e] != INVALID) {
1385 1385
          e = _right[e];
1386 1386
          while (_right[e] != INVALID) {
1387 1387
            e = _right[e];
1388 1388
          }
1389 1389
          Arc s = _parent[e];
1390
          _right.set(_parent[e], _left[e]);
1390
          _right[_parent[e]] = _left[e];
1391 1391
          if (_left[e] != INVALID) {
1392
            _parent.set(_left[e], _parent[e]);
1392
            _parent[_left[e]] = _parent[e];
1393 1393
          }
1394 1394

	
1395
          _left.set(e, _left[arc]);
1396
          _parent.set(_left[arc], e);
1397
          _right.set(e, _right[arc]);
1398
          _parent.set(_right[arc], e);
1395
          _left[e] = _left[arc];
1396
          _parent[_left[arc]] = e;
1397
          _right[e] = _right[arc];
1398
          _parent[_right[arc]] = e;
1399 1399

	
1400
          _parent.set(e, _parent[arc]);
1400
          _parent[e] = _parent[arc];
1401 1401
          if (_parent[arc] != INVALID) {
1402 1402
            if (_left[_parent[arc]] == arc) {
1403
              _left.set(_parent[arc], e);
1403
              _left[_parent[arc]] = e;
1404 1404
            } else {
1405
              _right.set(_parent[arc], e);
1405
              _right[_parent[arc]] = e;
1406 1406
            }
1407 1407
          }
1408 1408
          splay(s);
1409 1409
        } else {
1410
          _right.set(e, _right[arc]);
1411
          _parent.set(_right[arc], e);
1412
          _parent.set(e, _parent[arc]);
1410
          _right[e] = _right[arc];
1411
          _parent[_right[arc]] = e;
1412
          _parent[e] = _parent[arc];
1413 1413

	
1414 1414
          if (_parent[arc] != INVALID) {
1415 1415
            if (_left[_parent[arc]] == arc) {
1416
              _left.set(_parent[arc], e);
1416
              _left[_parent[arc]] = e;
1417 1417
            } else {
1418
              _right.set(_parent[arc], e);
1418
              _right[_parent[arc]] = e;
1419 1419
            }
1420 1420
          } else {
1421
            _head.set(_g.source(arc), e);
1421
            _head[_g.source(arc)] = e;
1422 1422
          }
1423 1423
        }
1424 1424
      }
1425 1425
    }
1426 1426

	
1427 1427
    Arc refreshRec(std::vector<Arc> &v,int a,int b)
1428 1428
    {
1429 1429
      int m=(a+b)/2;
1430 1430
      Arc me=v[m];
1431 1431
      if (a < m) {
1432 1432
        Arc left = refreshRec(v,a,m-1);
1433
        _left.set(me, left);
1434
        _parent.set(left, me);
1433
        _left[me] = left;
1434
        _parent[left] = me;
1435 1435
      } else {
1436
        _left.set(me, INVALID);
1436
        _left[me] = INVALID;
1437 1437
      }
1438 1438
      if (m < b) {
1439 1439
        Arc right = refreshRec(v,m+1,b);
1440
        _right.set(me, right);
1441
        _parent.set(right, me);
1440
        _right[me] = right;
1441
        _parent[right] = me;
1442 1442
      } else {
1443
        _right.set(me, INVALID);
1443
        _right[me] = INVALID;
1444 1444
      }
1445 1445
      return me;
1446 1446
    }
1447 1447

	
1448 1448
    void refresh() {
1449 1449
      for(NodeIt n(_g);n!=INVALID;++n) {
1450 1450
        std::vector<Arc> v;
1451 1451
        for(OutArcIt a(_g,n);a!=INVALID;++a) v.push_back(a);
1452 1452
        if (!v.empty()) {
1453 1453
          std::sort(v.begin(),v.end(),ArcLess(_g));
1454 1454
          Arc head = refreshRec(v,0,v.size()-1);
1455
          _head.set(n, head);
1456
          _parent.set(head, INVALID);
1455
          _head[n] = head;
1456
          _parent[head] = INVALID;
1457 1457
        }
1458
        else _head.set(n, INVALID);
1458
        else _head[n] = INVALID;
1459 1459
      }
1460 1460
    }
1461 1461

	
1462 1462
    void zig(Arc v) {
1463 1463
      Arc w = _parent[v];
1464
      _parent.set(v, _parent[w]);
1465
      _parent.set(w, v);
1466
      _left.set(w, _right[v]);
1467
      _right.set(v, w);
1464
      _parent[v] = _parent[w];
1465
      _parent[w] = v;
1466
      _left[w] = _right[v];
1467
      _right[v] = w;
1468 1468
      if (_parent[v] != INVALID) {
1469 1469
        if (_right[_parent[v]] == w) {
1470
          _right.set(_parent[v], v);
1470
          _right[_parent[v]] = v;
1471 1471
        } else {
1472
          _left.set(_parent[v], v);
1472
          _left[_parent[v]] = v;
1473 1473
        }
1474 1474
      }
1475 1475
      if (_left[w] != INVALID){
1476
        _parent.set(_left[w], w);
1476
        _parent[_left[w]] = w;
1477 1477
      }
1478 1478
    }
1479 1479

	
1480 1480
    void zag(Arc v) {
1481 1481
      Arc w = _parent[v];
1482
      _parent.set(v, _parent[w]);
1483
      _parent.set(w, v);
1484
      _right.set(w, _left[v]);
1485
      _left.set(v, w);
1482
      _parent[v] = _parent[w];
1483
      _parent[w] = v;
1484
      _right[w] = _left[v];
1485
      _left[v] = w;
1486 1486
      if (_parent[v] != INVALID){
1487 1487
        if (_left[_parent[v]] == w) {
1488
          _left.set(_parent[v], v);
1488
          _left[_parent[v]] = v;
1489 1489
        } else {
1490
          _right.set(_parent[v], v);
1490
          _right[_parent[v]] = v;
1491 1491
        }
1492 1492
      }
1493 1493
      if (_right[w] != INVALID){
1494
        _parent.set(_right[w], w);
1494
        _parent[_right[w]] = w;
1495 1495
      }
1496 1496
    }
1497 1497

	
1498 1498
    void splay(Arc v) {
1499 1499
      while (_parent[v] != INVALID) {
1500 1500
        if (v == _left[_parent[v]]) {
1501 1501
          if (_parent[_parent[v]] == INVALID) {
1502 1502
            zig(v);
1503 1503
          } else {
1504 1504
            if (_parent[v] == _left[_parent[_parent[v]]]) {
1505 1505
              zig(_parent[v]);
1506 1506
              zig(v);
1507 1507
            } else {
1508 1508
              zig(v);
1509 1509
              zag(v);
1510 1510
            }
1511 1511
          }
1512 1512
        } else {
1513 1513
          if (_parent[_parent[v]] == INVALID) {
1514 1514
            zag(v);
1515 1515
          } else {
1516 1516
            if (_parent[v] == _left[_parent[_parent[v]]]) {
1517 1517
              zag(v);
1518 1518
              zig(v);
1519 1519
            } else {
1520 1520
              zag(_parent[v]);
1521 1521
              zag(v);
1522 1522
            }
1523 1523
          }
1524 1524
        }
1525 1525
      }
1526 1526
      _head[_g.source(v)] = v;
1527 1527
    }
1528 1528

	
1529 1529

	
1530 1530
  public:
1531 1531

	
1532 1532
    ///Find an arc between two nodes.
1533 1533

	
1534 1534
    ///Find an arc between two nodes.
1535 1535
    ///\param s The source node.
1536 1536
    ///\param t The target node.
1537 1537
    ///\param p The previous arc between \c s and \c t. It it is INVALID or
1538 1538
    ///not given, the operator finds the first appropriate arc.
1539 1539
    ///\return An arc from \c s to \c t after \c p or
1540 1540
    ///\ref INVALID if there is no more.
1541 1541
    ///
1542 1542
    ///For example, you can count the number of arcs from \c u to \c v in the
Show white space 96 line context
... ...
@@ -31,113 +31,113 @@
31 31
#include <lemon/bits/traits.h>
32 32

	
33 33
namespace lemon {
34 34

	
35 35
  ///Class for handling "labels" in push-relabel type algorithms.
36 36

	
37 37
  ///A class for handling "labels" in push-relabel type algorithms.
38 38
  ///
39 39
  ///\ingroup auxdat
40 40
  ///Using this class you can assign "labels" (nonnegative integer numbers)
41 41
  ///to the edges or nodes of a graph, manipulate and query them through
42 42
  ///operations typically arising in "push-relabel" type algorithms.
43 43
  ///
44 44
  ///Each item is either \em active or not, and you can also choose a
45 45
  ///highest level active item.
46 46
  ///
47 47
  ///\sa LinkedElevator
48 48
  ///
49 49
  ///\param GR Type of the underlying graph.
50 50
  ///\param Item Type of the items the data is assigned to (\c GR::Node,
51 51
  ///\c GR::Arc or \c GR::Edge).
52 52
  template<class GR, class Item>
53 53
  class Elevator
54 54
  {
55 55
  public:
56 56

	
57 57
    typedef Item Key;
58 58
    typedef int Value;
59 59

	
60 60
  private:
61 61

	
62 62
    typedef Item *Vit;
63 63
    typedef typename ItemSetTraits<GR,Item>::template Map<Vit>::Type VitMap;
64 64
    typedef typename ItemSetTraits<GR,Item>::template Map<int>::Type IntMap;
65 65

	
66 66
    const GR &_g;
67 67
    int _max_level;
68 68
    int _item_num;
69 69
    VitMap _where;
70 70
    IntMap _level;
71 71
    std::vector<Item> _items;
72 72
    std::vector<Vit> _first;
73 73
    std::vector<Vit> _last_active;
74 74

	
75 75
    int _highest_active;
76 76

	
77 77
    void copy(Item i, Vit p)
78 78
    {
79
      _where.set(*p=i,p);
79
      _where[*p=i] = p;
80 80
    }
81 81
    void copy(Vit s, Vit p)
82 82
    {
83 83
      if(s!=p)
84 84
        {
85 85
          Item i=*s;
86 86
          *p=i;
87
          _where.set(i,p);
87
          _where[i] = p;
88 88
        }
89 89
    }
90 90
    void swap(Vit i, Vit j)
91 91
    {
92 92
      Item ti=*i;
93 93
      Vit ct = _where[ti];
94
      _where.set(ti,_where[*i=*j]);
95
      _where.set(*j,ct);
94
      _where[ti] = _where[*i=*j];
95
      _where[*j] = ct;
96 96
      *j=ti;
97 97
    }
98 98

	
99 99
  public:
100 100

	
101 101
    ///Constructor with given maximum level.
102 102

	
103 103
    ///Constructor with given maximum level.
104 104
    ///
105 105
    ///\param graph The underlying graph.
106 106
    ///\param max_level The maximum allowed level.
107 107
    ///Set the range of the possible labels to <tt>[0..max_level]</tt>.
108 108
    Elevator(const GR &graph,int max_level) :
109 109
      _g(graph),
110 110
      _max_level(max_level),
111 111
      _item_num(_max_level),
112 112
      _where(graph),
113 113
      _level(graph,0),
114 114
      _items(_max_level),
115 115
      _first(_max_level+2),
116 116
      _last_active(_max_level+2),
117 117
      _highest_active(-1) {}
118 118
    ///Constructor.
119 119

	
120 120
    ///Constructor.
121 121
    ///
122 122
    ///\param graph The underlying graph.
123 123
    ///Set the range of the possible labels to <tt>[0..max_level]</tt>,
124 124
    ///where \c max_level is equal to the number of labeled items in the graph.
125 125
    Elevator(const GR &graph) :
126 126
      _g(graph),
127 127
      _max_level(countItems<GR, Item>(graph)),
128 128
      _item_num(_max_level),
129 129
      _where(graph),
130 130
      _level(graph,0),
131 131
      _items(_max_level),
132 132
      _first(_max_level+2),
133 133
      _last_active(_max_level+2),
134 134
      _highest_active(-1)
135 135
    {
136 136
    }
137 137

	
138 138
    ///Activate item \c i.
139 139

	
140 140
    ///Activate item \c i.
141 141
    ///\pre Item \c i shouldn't be active before.
142 142
    void activate(Item i)
143 143
    {
... ...
@@ -181,314 +181,314 @@
181 181
    }
182 182
    ///Return the number of active items on level \c l.
183 183
    int activesOnLevel(int l) const
184 184
    {
185 185
      return _last_active[l]-_first[l]+1;
186 186
    }
187 187
    ///Return true if there is no active item on level \c l.
188 188
    bool activeFree(int l) const
189 189
    {
190 190
      return _last_active[l]<_first[l];
191 191
    }
192 192
    ///Return the maximum allowed level.
193 193
    int maxLevel() const
194 194
    {
195 195
      return _max_level;
196 196
    }
197 197

	
198 198
    ///\name Highest Active Item
199 199
    ///Functions for working with the highest level
200 200
    ///active item.
201 201

	
202 202
    ///@{
203 203

	
204 204
    ///Return a highest level active item.
205 205

	
206 206
    ///Return a highest level active item or INVALID if there is no active
207 207
    ///item.
208 208
    Item highestActive() const
209 209
    {
210 210
      return _highest_active>=0?*_last_active[_highest_active]:INVALID;
211 211
    }
212 212

	
213 213
    ///Return the highest active level.
214 214

	
215 215
    ///Return the level of the highest active item or -1 if there is no active
216 216
    ///item.
217 217
    int highestActiveLevel() const
218 218
    {
219 219
      return _highest_active;
220 220
    }
221 221

	
222 222
    ///Lift the highest active item by one.
223 223

	
224 224
    ///Lift the item returned by highestActive() by one.
225 225
    ///
226 226
    void liftHighestActive()
227 227
    {
228 228
      Item it = *_last_active[_highest_active];
229
      _level.set(it,_level[it]+1);
229
      ++_level[it];
230 230
      swap(_last_active[_highest_active]--,_last_active[_highest_active+1]);
231 231
      --_first[++_highest_active];
232 232
    }
233 233

	
234 234
    ///Lift the highest active item to the given level.
235 235

	
236 236
    ///Lift the item returned by highestActive() to level \c new_level.
237 237
    ///
238 238
    ///\warning \c new_level must be strictly higher
239 239
    ///than the current level.
240 240
    ///
241 241
    void liftHighestActive(int new_level)
242 242
    {
243 243
      const Item li = *_last_active[_highest_active];
244 244

	
245 245
      copy(--_first[_highest_active+1],_last_active[_highest_active]--);
246 246
      for(int l=_highest_active+1;l<new_level;l++)
247 247
        {
248 248
          copy(--_first[l+1],_first[l]);
249 249
          --_last_active[l];
250 250
        }
251 251
      copy(li,_first[new_level]);
252
      _level.set(li,new_level);
252
      _level[li] = new_level;
253 253
      _highest_active=new_level;
254 254
    }
255 255

	
256 256
    ///Lift the highest active item to the top level.
257 257

	
258 258
    ///Lift the item returned by highestActive() to the top level and
259 259
    ///deactivate it.
260 260
    void liftHighestActiveToTop()
261 261
    {
262 262
      const Item li = *_last_active[_highest_active];
263 263

	
264 264
      copy(--_first[_highest_active+1],_last_active[_highest_active]--);
265 265
      for(int l=_highest_active+1;l<_max_level;l++)
266 266
        {
267 267
          copy(--_first[l+1],_first[l]);
268 268
          --_last_active[l];
269 269
        }
270 270
      copy(li,_first[_max_level]);
271 271
      --_last_active[_max_level];
272
      _level.set(li,_max_level);
272
      _level[li] = _max_level;
273 273

	
274 274
      while(_highest_active>=0 &&
275 275
            _last_active[_highest_active]<_first[_highest_active])
276 276
        _highest_active--;
277 277
    }
278 278

	
279 279
    ///@}
280 280

	
281 281
    ///\name Active Item on Certain Level
282 282
    ///Functions for working with the active items.
283 283

	
284 284
    ///@{
285 285

	
286 286
    ///Return an active item on level \c l.
287 287

	
288 288
    ///Return an active item on level \c l or \ref INVALID if there is no such
289 289
    ///an item. (\c l must be from the range [0...\c max_level].
290 290
    Item activeOn(int l) const
291 291
    {
292 292
      return _last_active[l]>=_first[l]?*_last_active[l]:INVALID;
293 293
    }
294 294

	
295 295
    ///Lift the active item returned by \c activeOn(level) by one.
296 296

	
297 297
    ///Lift the active item returned by \ref activeOn() "activeOn(level)"
298 298
    ///by one.
299 299
    Item liftActiveOn(int level)
300 300
    {
301 301
      Item it =*_last_active[level];
302
      _level.set(it,_level[it]+1);
302
      ++_level[it];
303 303
      swap(_last_active[level]--, --_first[level+1]);
304 304
      if (level+1>_highest_active) ++_highest_active;
305 305
    }
306 306

	
307 307
    ///Lift the active item returned by \c activeOn(level) to the given level.
308 308

	
309 309
    ///Lift the active item returned by \ref activeOn() "activeOn(level)"
310 310
    ///to the given level.
311 311
    void liftActiveOn(int level, int new_level)
312 312
    {
313 313
      const Item ai = *_last_active[level];
314 314

	
315 315
      copy(--_first[level+1], _last_active[level]--);
316 316
      for(int l=level+1;l<new_level;l++)
317 317
        {
318 318
          copy(_last_active[l],_first[l]);
319 319
          copy(--_first[l+1], _last_active[l]--);
320 320
        }
321 321
      copy(ai,_first[new_level]);
322
      _level.set(ai,new_level);
322
      _level[ai] = new_level;
323 323
      if (new_level>_highest_active) _highest_active=new_level;
324 324
    }
325 325

	
326 326
    ///Lift the active item returned by \c activeOn(level) to the top level.
327 327

	
328 328
    ///Lift the active item returned by \ref activeOn() "activeOn(level)"
329 329
    ///to the top level and deactivate it.
330 330
    void liftActiveToTop(int level)
331 331
    {
332 332
      const Item ai = *_last_active[level];
333 333

	
334 334
      copy(--_first[level+1],_last_active[level]--);
335 335
      for(int l=level+1;l<_max_level;l++)
336 336
        {
337 337
          copy(_last_active[l],_first[l]);
338 338
          copy(--_first[l+1], _last_active[l]--);
339 339
        }
340 340
      copy(ai,_first[_max_level]);
341 341
      --_last_active[_max_level];
342
      _level.set(ai,_max_level);
342
      _level[ai] = _max_level;
343 343

	
344 344
      if (_highest_active==level) {
345 345
        while(_highest_active>=0 &&
346 346
              _last_active[_highest_active]<_first[_highest_active])
347 347
          _highest_active--;
348 348
      }
349 349
    }
350 350

	
351 351
    ///@}
352 352

	
353 353
    ///Lift an active item to a higher level.
354 354

	
355 355
    ///Lift an active item to a higher level.
356 356
    ///\param i The item to be lifted. It must be active.
357 357
    ///\param new_level The new level of \c i. It must be strictly higher
358 358
    ///than the current level.
359 359
    ///
360 360
    void lift(Item i, int new_level)
361 361
    {
362 362
      const int lo = _level[i];
363 363
      const Vit w = _where[i];
364 364

	
365 365
      copy(_last_active[lo],w);
366 366
      copy(--_first[lo+1],_last_active[lo]--);
367 367
      for(int l=lo+1;l<new_level;l++)
368 368
        {
369 369
          copy(_last_active[l],_first[l]);
370 370
          copy(--_first[l+1],_last_active[l]--);
371 371
        }
372 372
      copy(i,_first[new_level]);
373
      _level.set(i,new_level);
373
      _level[i] = new_level;
374 374
      if(new_level>_highest_active) _highest_active=new_level;
375 375
    }
376 376

	
377 377
    ///Move an inactive item to the top but one level (in a dirty way).
378 378

	
379 379
    ///This function moves an inactive item from the top level to the top
380 380
    ///but one level (in a dirty way).
381 381
    ///\warning It makes the underlying datastructure corrupt, so use it
382 382
    ///only if you really know what it is for.
383 383
    ///\pre The item is on the top level.
384 384
    void dirtyTopButOne(Item i) {
385
      _level.set(i,_max_level - 1);
385
      _level[i] = _max_level - 1;
386 386
    }
387 387

	
388 388
    ///Lift all items on and above the given level to the top level.
389 389

	
390 390
    ///This function lifts all items on and above level \c l to the top
391 391
    ///level and deactivates them.
392 392
    void liftToTop(int l)
393 393
    {
394 394
      const Vit f=_first[l];
395 395
      const Vit tl=_first[_max_level];
396 396
      for(Vit i=f;i!=tl;++i)
397
        _level.set(*i,_max_level);
397
        _level[*i] = _max_level;
398 398
      for(int i=l;i<=_max_level;i++)
399 399
        {
400 400
          _first[i]=f;
401 401
          _last_active[i]=f-1;
402 402
        }
403 403
      for(_highest_active=l-1;
404 404
          _highest_active>=0 &&
405 405
            _last_active[_highest_active]<_first[_highest_active];
406 406
          _highest_active--) ;
407 407
    }
408 408

	
409 409
  private:
410 410
    int _init_lev;
411 411
    Vit _init_num;
412 412

	
413 413
  public:
414 414

	
415 415
    ///\name Initialization
416 416
    ///Using these functions you can initialize the levels of the items.
417 417
    ///\n
418 418
    ///The initialization must be started with calling \c initStart().
419 419
    ///Then the items should be listed level by level starting with the
420 420
    ///lowest one (level 0) using \c initAddItem() and \c initNewLevel().
421 421
    ///Finally \c initFinish() must be called.
422 422
    ///The items not listed are put on the highest level.
423 423
    ///@{
424 424

	
425 425
    ///Start the initialization process.
426 426
    void initStart()
427 427
    {
428 428
      _init_lev=0;
429 429
      _init_num=&_items[0];
430 430
      _first[0]=&_items[0];
431 431
      _last_active[0]=&_items[0]-1;
432 432
      Vit n=&_items[0];
433 433
      for(typename ItemSetTraits<GR,Item>::ItemIt i(_g);i!=INVALID;++i)
434 434
        {
435 435
          *n=i;
436
          _where.set(i,n);
437
          _level.set(i,_max_level);
436
          _where[i] = n;
437
          _level[i] = _max_level;
438 438
          ++n;
439 439
        }
440 440
    }
441 441

	
442 442
    ///Add an item to the current level.
443 443
    void initAddItem(Item i)
444 444
    {
445 445
      swap(_where[i],_init_num);
446
      _level.set(i,_init_lev);
446
      _level[i] = _init_lev;
447 447
      ++_init_num;
448 448
    }
449 449

	
450 450
    ///Start a new level.
451 451

	
452 452
    ///Start a new level.
453 453
    ///It shouldn't be used before the items on level 0 are listed.
454 454
    void initNewLevel()
455 455
    {
456 456
      _init_lev++;
457 457
      _first[_init_lev]=_init_num;
458 458
      _last_active[_init_lev]=_init_num-1;
459 459
    }
460 460

	
461 461
    ///Finalize the initialization process.
462 462
    void initFinish()
463 463
    {
464 464
      for(_init_lev++;_init_lev<=_max_level;_init_lev++)
465 465
        {
466 466
          _first[_init_lev]=_init_num;
467 467
          _last_active[_init_lev]=_init_num-1;
468 468
        }
469 469
      _first[_max_level+1]=&_items[0]+_item_num;
470 470
      _last_active[_max_level+1]=&_items[0]+_item_num-1;
471 471
      _highest_active = -1;
472 472
    }
473 473

	
474 474
    ///@}
475 475

	
476 476
  };
477 477

	
478 478
  ///Class for handling "labels" in push-relabel type algorithms.
479 479

	
480 480
  ///A class for handling "labels" in push-relabel type algorithms.
481 481
  ///
482 482
  ///\ingroup auxdat
483 483
  ///Using this class you can assign "labels" (nonnegative integer numbers)
484 484
  ///to the edges or nodes of a graph, manipulate and query them through
485 485
  ///operations typically arising in "push-relabel" type algorithms.
486 486
  ///
487 487
  ///Each item is either \em active or not, and you can also choose a
488 488
  ///highest level active item.
489 489
  ///
490 490
  ///\sa Elevator
491 491
  ///
492 492
  ///\param GR Type of the underlying graph.
493 493
  ///\param Item Type of the items the data is assigned to (\c GR::Node,
494 494
  ///\c GR::Arc or \c GR::Edge).
... ...
@@ -506,477 +506,477 @@
506 506
    typedef typename ItemSetTraits<GR,Item>::
507 507
    template Map<int>::Type IntMap;
508 508
    typedef typename ItemSetTraits<GR,Item>::
509 509
    template Map<bool>::Type BoolMap;
510 510

	
511 511
    const GR &_graph;
512 512
    int _max_level;
513 513
    int _item_num;
514 514
    std::vector<Item> _first, _last;
515 515
    ItemMap _prev, _next;
516 516
    int _highest_active;
517 517
    IntMap _level;
518 518
    BoolMap _active;
519 519

	
520 520
  public:
521 521
    ///Constructor with given maximum level.
522 522

	
523 523
    ///Constructor with given maximum level.
524 524
    ///
525 525
    ///\param graph The underlying graph.
526 526
    ///\param max_level The maximum allowed level.
527 527
    ///Set the range of the possible labels to <tt>[0..max_level]</tt>.
528 528
    LinkedElevator(const GR& graph, int max_level)
529 529
      : _graph(graph), _max_level(max_level), _item_num(_max_level),
530 530
        _first(_max_level + 1), _last(_max_level + 1),
531 531
        _prev(graph), _next(graph),
532 532
        _highest_active(-1), _level(graph), _active(graph) {}
533 533

	
534 534
    ///Constructor.
535 535

	
536 536
    ///Constructor.
537 537
    ///
538 538
    ///\param graph The underlying graph.
539 539
    ///Set the range of the possible labels to <tt>[0..max_level]</tt>,
540 540
    ///where \c max_level is equal to the number of labeled items in the graph.
541 541
    LinkedElevator(const GR& graph)
542 542
      : _graph(graph), _max_level(countItems<GR, Item>(graph)),
543 543
        _item_num(_max_level),
544 544
        _first(_max_level + 1), _last(_max_level + 1),
545 545
        _prev(graph, INVALID), _next(graph, INVALID),
546 546
        _highest_active(-1), _level(graph), _active(graph) {}
547 547

	
548 548

	
549 549
    ///Activate item \c i.
550 550

	
551 551
    ///Activate item \c i.
552 552
    ///\pre Item \c i shouldn't be active before.
553 553
    void activate(Item i) {
554
      _active.set(i, true);
554
      _active[i] = true;
555 555

	
556 556
      int level = _level[i];
557 557
      if (level > _highest_active) {
558 558
        _highest_active = level;
559 559
      }
560 560

	
561 561
      if (_prev[i] == INVALID || _active[_prev[i]]) return;
562 562
      //unlace
563
      _next.set(_prev[i], _next[i]);
563
      _next[_prev[i]] = _next[i];
564 564
      if (_next[i] != INVALID) {
565
        _prev.set(_next[i], _prev[i]);
565
        _prev[_next[i]] = _prev[i];
566 566
      } else {
567 567
        _last[level] = _prev[i];
568 568
      }
569 569
      //lace
570
      _next.set(i, _first[level]);
571
      _prev.set(_first[level], i);
572
      _prev.set(i, INVALID);
570
      _next[i] = _first[level];
571
      _prev[_first[level]] = i;
572
      _prev[i] = INVALID;
573 573
      _first[level] = i;
574 574

	
575 575
    }
576 576

	
577 577
    ///Deactivate item \c i.
578 578

	
579 579
    ///Deactivate item \c i.
580 580
    ///\pre Item \c i must be active before.
581 581
    void deactivate(Item i) {
582
      _active.set(i, false);
582
      _active[i] = false;
583 583
      int level = _level[i];
584 584

	
585 585
      if (_next[i] == INVALID || !_active[_next[i]])
586 586
        goto find_highest_level;
587 587

	
588 588
      //unlace
589
      _prev.set(_next[i], _prev[i]);
589
      _prev[_next[i]] = _prev[i];
590 590
      if (_prev[i] != INVALID) {
591
        _next.set(_prev[i], _next[i]);
591
        _next[_prev[i]] = _next[i];
592 592
      } else {
593 593
        _first[_level[i]] = _next[i];
594 594
      }
595 595
      //lace
596
      _prev.set(i, _last[level]);
597
      _next.set(_last[level], i);
598
      _next.set(i, INVALID);
596
      _prev[i] = _last[level];
597
      _next[_last[level]] = i;
598
      _next[i] = INVALID;
599 599
      _last[level] = i;
600 600

	
601 601
    find_highest_level:
602 602
      if (level == _highest_active) {
603 603
        while (_highest_active >= 0 && activeFree(_highest_active))
604 604
          --_highest_active;
605 605
      }
606 606
    }
607 607

	
608 608
    ///Query whether item \c i is active
609 609
    bool active(Item i) const { return _active[i]; }
610 610

	
611 611
    ///Return the level of item \c i.
612 612
    int operator[](Item i) const { return _level[i]; }
613 613

	
614 614
    ///Return the number of items on level \c l.
615 615
    int onLevel(int l) const {
616 616
      int num = 0;
617 617
      Item n = _first[l];
618 618
      while (n != INVALID) {
619 619
        ++num;
620 620
        n = _next[n];
621 621
      }
622 622
      return num;
623 623
    }
624 624

	
625 625
    ///Return true if the level is empty.
626 626
    bool emptyLevel(int l) const {
627 627
      return _first[l] == INVALID;
628 628
    }
629 629

	
630 630
    ///Return the number of items above level \c l.
631 631
    int aboveLevel(int l) const {
632 632
      int num = 0;
633 633
      for (int level = l + 1; level < _max_level; ++level)
634 634
        num += onLevel(level);
635 635
      return num;
636 636
    }
637 637

	
638 638
    ///Return the number of active items on level \c l.
639 639
    int activesOnLevel(int l) const {
640 640
      int num = 0;
641 641
      Item n = _first[l];
642 642
      while (n != INVALID && _active[n]) {
643 643
        ++num;
644 644
        n = _next[n];
645 645
      }
646 646
      return num;
647 647
    }
648 648

	
649 649
    ///Return true if there is no active item on level \c l.
650 650
    bool activeFree(int l) const {
651 651
      return _first[l] == INVALID || !_active[_first[l]];
652 652
    }
653 653

	
654 654
    ///Return the maximum allowed level.
655 655
    int maxLevel() const {
656 656
      return _max_level;
657 657
    }
658 658

	
659 659
    ///\name Highest Active Item
660 660
    ///Functions for working with the highest level
661 661
    ///active item.
662 662

	
663 663
    ///@{
664 664

	
665 665
    ///Return a highest level active item.
666 666

	
667 667
    ///Return a highest level active item or INVALID if there is no active
668 668
    ///item.
669 669
    Item highestActive() const {
670 670
      return _highest_active >= 0 ? _first[_highest_active] : INVALID;
671 671
    }
672 672

	
673 673
    ///Return the highest active level.
674 674

	
675 675
    ///Return the level of the highest active item or -1 if there is no active
676 676
    ///item.
677 677
    int highestActiveLevel() const {
678 678
      return _highest_active;
679 679
    }
680 680

	
681 681
    ///Lift the highest active item by one.
682 682

	
683 683
    ///Lift the item returned by highestActive() by one.
684 684
    ///
685 685
    void liftHighestActive() {
686 686
      Item i = _first[_highest_active];
687 687
      if (_next[i] != INVALID) {
688
        _prev.set(_next[i], INVALID);
688
        _prev[_next[i]] = INVALID;
689 689
        _first[_highest_active] = _next[i];
690 690
      } else {
691 691
        _first[_highest_active] = INVALID;
692 692
        _last[_highest_active] = INVALID;
693 693
      }
694
      _level.set(i, ++_highest_active);
694
      _level[i] = ++_highest_active;
695 695
      if (_first[_highest_active] == INVALID) {
696 696
        _first[_highest_active] = i;
697 697
        _last[_highest_active] = i;
698
        _prev.set(i, INVALID);
699
        _next.set(i, INVALID);
698
        _prev[i] = INVALID;
699
        _next[i] = INVALID;
700 700
      } else {
701
        _prev.set(_first[_highest_active], i);
702
        _next.set(i, _first[_highest_active]);
701
        _prev[_first[_highest_active]] = i;
702
        _next[i] = _first[_highest_active];
703 703
        _first[_highest_active] = i;
704 704
      }
705 705
    }
706 706

	
707 707
    ///Lift the highest active item to the given level.
708 708

	
709 709
    ///Lift the item returned by highestActive() to level \c new_level.
710 710
    ///
711 711
    ///\warning \c new_level must be strictly higher
712 712
    ///than the current level.
713 713
    ///
714 714
    void liftHighestActive(int new_level) {
715 715
      Item i = _first[_highest_active];
716 716
      if (_next[i] != INVALID) {
717
        _prev.set(_next[i], INVALID);
717
        _prev[_next[i]] = INVALID;
718 718
        _first[_highest_active] = _next[i];
719 719
      } else {
720 720
        _first[_highest_active] = INVALID;
721 721
        _last[_highest_active] = INVALID;
722 722
      }
723
      _level.set(i, _highest_active = new_level);
723
      _level[i] = _highest_active = new_level;
724 724
      if (_first[_highest_active] == INVALID) {
725 725
        _first[_highest_active] = _last[_highest_active] = i;
726
        _prev.set(i, INVALID);
727
        _next.set(i, INVALID);
726
        _prev[i] = INVALID;
727
        _next[i] = INVALID;
728 728
      } else {
729
        _prev.set(_first[_highest_active], i);
730
        _next.set(i, _first[_highest_active]);
729
        _prev[_first[_highest_active]] = i;
730
        _next[i] = _first[_highest_active];
731 731
        _first[_highest_active] = i;
732 732
      }
733 733
    }
734 734

	
735 735
    ///Lift the highest active item to the top level.
736 736

	
737 737
    ///Lift the item returned by highestActive() to the top level and
738 738
    ///deactivate it.
739 739
    void liftHighestActiveToTop() {
740 740
      Item i = _first[_highest_active];
741
      _level.set(i, _max_level);
741
      _level[i] = _max_level;
742 742
      if (_next[i] != INVALID) {
743
        _prev.set(_next[i], INVALID);
743
        _prev[_next[i]] = INVALID;
744 744
        _first[_highest_active] = _next[i];
745 745
      } else {
746 746
        _first[_highest_active] = INVALID;
747 747
        _last[_highest_active] = INVALID;
748 748
      }
749 749
      while (_highest_active >= 0 && activeFree(_highest_active))
750 750
        --_highest_active;
751 751
    }
752 752

	
753 753
    ///@}
754 754

	
755 755
    ///\name Active Item on Certain Level
756 756
    ///Functions for working with the active items.
757 757

	
758 758
    ///@{
759 759

	
760 760
    ///Return an active item on level \c l.
761 761

	
762 762
    ///Return an active item on level \c l or \ref INVALID if there is no such
763 763
    ///an item. (\c l must be from the range [0...\c max_level].
764 764
    Item activeOn(int l) const
765 765
    {
766 766
      return _active[_first[l]] ? _first[l] : INVALID;
767 767
    }
768 768

	
769 769
    ///Lift the active item returned by \c activeOn(l) by one.
770 770

	
771 771
    ///Lift the active item returned by \ref activeOn() "activeOn(l)"
772 772
    ///by one.
773 773
    Item liftActiveOn(int l)
774 774
    {
775 775
      Item i = _first[l];
776 776
      if (_next[i] != INVALID) {
777
        _prev.set(_next[i], INVALID);
777
        _prev[_next[i]] = INVALID;
778 778
        _first[l] = _next[i];
779 779
      } else {
780 780
        _first[l] = INVALID;
781 781
        _last[l] = INVALID;
782 782
      }
783
      _level.set(i, ++l);
783
      _level[i] = ++l;
784 784
      if (_first[l] == INVALID) {
785 785
        _first[l] = _last[l] = i;
786
        _prev.set(i, INVALID);
787
        _next.set(i, INVALID);
786
        _prev[i] = INVALID;
787
        _next[i] = INVALID;
788 788
      } else {
789
        _prev.set(_first[l], i);
790
        _next.set(i, _first[l]);
789
        _prev[_first[l]] = i;
790
        _next[i] = _first[l];
791 791
        _first[l] = i;
792 792
      }
793 793
      if (_highest_active < l) {
794 794
        _highest_active = l;
795 795
      }
796 796
    }
797 797

	
798 798
    ///Lift the active item returned by \c activeOn(l) to the given level.
799 799

	
800 800
    ///Lift the active item returned by \ref activeOn() "activeOn(l)"
801 801
    ///to the given level.
802 802
    void liftActiveOn(int l, int new_level)
803 803
    {
804 804
      Item i = _first[l];
805 805
      if (_next[i] != INVALID) {
806
        _prev.set(_next[i], INVALID);
806
        _prev[_next[i]] = INVALID;
807 807
        _first[l] = _next[i];
808 808
      } else {
809 809
        _first[l] = INVALID;
810 810
        _last[l] = INVALID;
811 811
      }
812
      _level.set(i, l = new_level);
812
      _level[i] = l = new_level;
813 813
      if (_first[l] == INVALID) {
814 814
        _first[l] = _last[l] = i;
815
        _prev.set(i, INVALID);
816
        _next.set(i, INVALID);
815
        _prev[i] = INVALID;
816
        _next[i] = INVALID;
817 817
      } else {
818
        _prev.set(_first[l], i);
819
        _next.set(i, _first[l]);
818
        _prev[_first[l]] = i;
819
        _next[i] = _first[l];
820 820
        _first[l] = i;
821 821
      }
822 822
      if (_highest_active < l) {
823 823
        _highest_active = l;
824 824
      }
825 825
    }
826 826

	
827 827
    ///Lift the active item returned by \c activeOn(l) to the top level.
828 828

	
829 829
    ///Lift the active item returned by \ref activeOn() "activeOn(l)"
830 830
    ///to the top level and deactivate it.
831 831
    void liftActiveToTop(int l)
832 832
    {
833 833
      Item i = _first[l];
834 834
      if (_next[i] != INVALID) {
835
        _prev.set(_next[i], INVALID);
835
        _prev[_next[i]] = INVALID;
836 836
        _first[l] = _next[i];
837 837
      } else {
838 838
        _first[l] = INVALID;
839 839
        _last[l] = INVALID;
840 840
      }
841
      _level.set(i, _max_level);
841
      _level[i] = _max_level;
842 842
      if (l == _highest_active) {
843 843
        while (_highest_active >= 0 && activeFree(_highest_active))
844 844
          --_highest_active;
845 845
      }
846 846
    }
847 847

	
848 848
    ///@}
849 849

	
850 850
    /// \brief Lift an active item to a higher level.
851 851
    ///
852 852
    /// Lift an active item to a higher level.
853 853
    /// \param i The item to be lifted. It must be active.
854 854
    /// \param new_level The new level of \c i. It must be strictly higher
855 855
    /// than the current level.
856 856
    ///
857 857
    void lift(Item i, int new_level) {
858 858
      if (_next[i] != INVALID) {
859
        _prev.set(_next[i], _prev[i]);
859
        _prev[_next[i]] = _prev[i];
860 860
      } else {
861 861
        _last[new_level] = _prev[i];
862 862
      }
863 863
      if (_prev[i] != INVALID) {
864
        _next.set(_prev[i], _next[i]);
864
        _next[_prev[i]] = _next[i];
865 865
      } else {
866 866
        _first[new_level] = _next[i];
867 867
      }
868
      _level.set(i, new_level);
868
      _level[i] = new_level;
869 869
      if (_first[new_level] == INVALID) {
870 870
        _first[new_level] = _last[new_level] = i;
871
        _prev.set(i, INVALID);
872
        _next.set(i, INVALID);
871
        _prev[i] = INVALID;
872
        _next[i] = INVALID;
873 873
      } else {
874
        _prev.set(_first[new_level], i);
875
        _next.set(i, _first[new_level]);
874
        _prev[_first[new_level]] = i;
875
        _next[i] = _first[new_level];
876 876
        _first[new_level] = i;
877 877
      }
878 878
      if (_highest_active < new_level) {
879 879
        _highest_active = new_level;
880 880
      }
881 881
    }
882 882

	
883 883
    ///Move an inactive item to the top but one level (in a dirty way).
884 884

	
885 885
    ///This function moves an inactive item from the top level to the top
886 886
    ///but one level (in a dirty way).
887 887
    ///\warning It makes the underlying datastructure corrupt, so use it
888 888
    ///only if you really know what it is for.
889 889
    ///\pre The item is on the top level.
890 890
    void dirtyTopButOne(Item i) {
891
      _level.set(i, _max_level - 1);
891
      _level[i] = _max_level - 1;
892 892
    }
893 893

	
894 894
    ///Lift all items on and above the given level to the top level.
895 895

	
896 896
    ///This function lifts all items on and above level \c l to the top
897 897
    ///level and deactivates them.
898 898
    void liftToTop(int l)  {
899 899
      for (int i = l + 1; _first[i] != INVALID; ++i) {
900 900
        Item n = _first[i];
901 901
        while (n != INVALID) {
902
          _level.set(n, _max_level);
902
          _level[n] = _max_level;
903 903
          n = _next[n];
904 904
        }
905 905
        _first[i] = INVALID;
906 906
        _last[i] = INVALID;
907 907
      }
908 908
      if (_highest_active > l - 1) {
909 909
        _highest_active = l - 1;
910 910
        while (_highest_active >= 0 && activeFree(_highest_active))
911 911
          --_highest_active;
912 912
      }
913 913
    }
914 914

	
915 915
  private:
916 916

	
917 917
    int _init_level;
918 918

	
919 919
  public:
920 920

	
921 921
    ///\name Initialization
922 922
    ///Using these functions you can initialize the levels of the items.
923 923
    ///\n
924 924
    ///The initialization must be started with calling \c initStart().
925 925
    ///Then the items should be listed level by level starting with the
926 926
    ///lowest one (level 0) using \c initAddItem() and \c initNewLevel().
927 927
    ///Finally \c initFinish() must be called.
928 928
    ///The items not listed are put on the highest level.
929 929
    ///@{
930 930

	
931 931
    ///Start the initialization process.
932 932
    void initStart() {
933 933

	
934 934
      for (int i = 0; i <= _max_level; ++i) {
935 935
        _first[i] = _last[i] = INVALID;
936 936
      }
937 937
      _init_level = 0;
938 938
      for(typename ItemSetTraits<GR,Item>::ItemIt i(_graph);
939 939
          i != INVALID; ++i) {
940
        _level.set(i, _max_level);
941
        _active.set(i, false);
940
        _level[i] = _max_level;
941
        _active[i] = false;
942 942
      }
943 943
    }
944 944

	
945 945
    ///Add an item to the current level.
946 946
    void initAddItem(Item i) {
947
      _level.set(i, _init_level);
947
      _level[i] = _init_level;
948 948
      if (_last[_init_level] == INVALID) {
949 949
        _first[_init_level] = i;
950 950
        _last[_init_level] = i;
951
        _prev.set(i, INVALID);
952
        _next.set(i, INVALID);
951
        _prev[i] = INVALID;
952
        _next[i] = INVALID;
953 953
      } else {
954
        _prev.set(i, _last[_init_level]);
955
        _next.set(i, INVALID);
956
        _next.set(_last[_init_level], i);
954
        _prev[i] = _last[_init_level];
955
        _next[i] = INVALID;
956
        _next[_last[_init_level]] = i;
957 957
        _last[_init_level] = i;
958 958
      }
959 959
    }
960 960

	
961 961
    ///Start a new level.
962 962

	
963 963
    ///Start a new level.
964 964
    ///It shouldn't be used before the items on level 0 are listed.
965 965
    void initNewLevel() {
966 966
      ++_init_level;
967 967
    }
968 968

	
969 969
    ///Finalize the initialization process.
970 970
    void initFinish() {
971 971
      _highest_active = -1;
972 972
    }
973 973

	
974 974
    ///@}
975 975

	
976 976
  };
977 977

	
978 978

	
979 979
} //END OF NAMESPACE LEMON
980 980

	
981 981
#endif
982 982

	
Show white space 96 line context
... ...
@@ -98,144 +98,144 @@
98 98
      }
99 99
      if (!_order) {
100 100
	_order = new typename Graph::template NodeMap<int>(_graph);
101 101
      }
102 102
    }
103 103

	
104 104
    void destroyStructures() {
105 105
      if (_pred) {
106 106
	delete _pred;
107 107
      }
108 108
      if (_weight) {
109 109
	delete _weight;
110 110
      }
111 111
      if (_order) {
112 112
	delete _order;
113 113
      }
114 114
    }
115 115
  
116 116
  public:
117 117

	
118 118
    /// \brief Constructor
119 119
    ///
120 120
    /// Constructor
121 121
    /// \param graph The undirected graph the algorithm runs on.
122 122
    /// \param capacity The edge capacity map.
123 123
    GomoryHu(const Graph& graph, const Capacity& capacity) 
124 124
      : _graph(graph), _capacity(capacity),
125 125
	_pred(0), _weight(0), _order(0) 
126 126
    {
127 127
      checkConcept<concepts::ReadMap<Edge, Value>, Capacity>();
128 128
    }
129 129

	
130 130

	
131 131
    /// \brief Destructor
132 132
    ///
133 133
    /// Destructor
134 134
    ~GomoryHu() {
135 135
      destroyStructures();
136 136
    }
137 137

	
138 138
  private:
139 139
  
140 140
    // Initialize the internal data structures
141 141
    void init() {
142 142
      createStructures();
143 143

	
144 144
      _root = NodeIt(_graph);
145 145
      for (NodeIt n(_graph); n != INVALID; ++n) {
146
	_pred->set(n, _root);
147
	_order->set(n, -1);
146
        (*_pred)[n] = _root;
147
        (*_order)[n] = -1;
148 148
      }
149
      _pred->set(_root, INVALID);
150
      _weight->set(_root, std::numeric_limits<Value>::max()); 
149
      (*_pred)[_root] = INVALID;
150
      (*_weight)[_root] = std::numeric_limits<Value>::max(); 
151 151
    }
152 152

	
153 153

	
154 154
    // Start the algorithm
155 155
    void start() {
156 156
      Preflow<Graph, Capacity> fa(_graph, _capacity, _root, INVALID);
157 157

	
158 158
      for (NodeIt n(_graph); n != INVALID; ++n) {
159 159
	if (n == _root) continue;
160 160

	
161 161
	Node pn = (*_pred)[n];
162 162
	fa.source(n);
163 163
	fa.target(pn);
164 164

	
165 165
	fa.runMinCut();
166 166

	
167
	_weight->set(n, fa.flowValue());
167
	(*_weight)[n] = fa.flowValue();
168 168

	
169 169
	for (NodeIt nn(_graph); nn != INVALID; ++nn) {
170 170
	  if (nn != n && fa.minCut(nn) && (*_pred)[nn] == pn) {
171
	    _pred->set(nn, n);
171
	    (*_pred)[nn] = n;
172 172
	  }
173 173
	}
174 174
	if ((*_pred)[pn] != INVALID && fa.minCut((*_pred)[pn])) {
175
	  _pred->set(n, (*_pred)[pn]);
176
	  _pred->set(pn, n);
177
	  _weight->set(n, (*_weight)[pn]);
178
	  _weight->set(pn, fa.flowValue());	
175
	  (*_pred)[n] = (*_pred)[pn];
176
	  (*_pred)[pn] = n;
177
	  (*_weight)[n] = (*_weight)[pn];
178
	  (*_weight)[pn] = fa.flowValue();
179 179
	}
180 180
      }
181 181

	
182
      _order->set(_root, 0);
182
      (*_order)[_root] = 0;
183 183
      int index = 1;
184 184

	
185 185
      for (NodeIt n(_graph); n != INVALID; ++n) {
186 186
	std::vector<Node> st;
187 187
	Node nn = n;
188 188
	while ((*_order)[nn] == -1) {
189 189
	  st.push_back(nn);
190 190
	  nn = (*_pred)[nn];
191 191
	}
192 192
	while (!st.empty()) {
193
	  _order->set(st.back(), index++);
193
	  (*_order)[st.back()] = index++;
194 194
	  st.pop_back();
195 195
	}
196 196
      }
197 197
    }
198 198

	
199 199
  public:
200 200

	
201 201
    ///\name Execution Control
202 202
 
203 203
    ///@{
204 204

	
205 205
    /// \brief Run the Gomory-Hu algorithm.
206 206
    ///
207 207
    /// This function runs the Gomory-Hu algorithm.
208 208
    void run() {
209 209
      init();
210 210
      start();
211 211
    }
212 212
    
213 213
    /// @}
214 214

	
215 215
    ///\name Query Functions
216 216
    ///The results of the algorithm can be obtained using these
217 217
    ///functions.\n
218 218
    ///\ref run() "run()" should be called before using them.\n
219 219
    ///See also \ref MinCutNodeIt and \ref MinCutEdgeIt.
220 220

	
221 221
    ///@{
222 222

	
223 223
    /// \brief Return the predecessor node in the Gomory-Hu tree.
224 224
    ///
225 225
    /// This function returns the predecessor node in the Gomory-Hu tree.
226 226
    /// If the node is
227 227
    /// the root of the Gomory-Hu tree, then it returns \c INVALID.
228 228
    Node predNode(const Node& node) {
229 229
      return (*_pred)[node];
230 230
    }
231 231

	
232 232
    /// \brief Return the distance from the root node in the Gomory-Hu tree.
233 233
    ///
234 234
    /// This function returns the distance of \c node from the root node
235 235
    /// in the Gomory-Hu tree.
236 236
    int rootDist(const Node& node) {
237 237
      return (*_order)[node];
238 238
    }
239 239

	
240 240
    /// \brief Return the weight of the predecessor edge in the
241 241
    /// Gomory-Hu tree.
... ...
@@ -264,99 +264,99 @@
264 264
	  if ((*_weight)[sn] <= value) value = (*_weight)[sn];
265 265
	  sn = (*_pred)[sn];
266 266
	}
267 267
      }
268 268
      return value;
269 269
    }
270 270

	
271 271
    /// \brief Return the minimum cut between two nodes
272 272
    ///
273 273
    /// This function returns the minimum cut between the nodes \c s and \c t
274 274
    /// in the \c cutMap parameter by setting the nodes in the component of
275 275
    /// \c s to \c true and the other nodes to \c false.
276 276
    ///
277 277
    /// For higher level interfaces, see MinCutNodeIt and MinCutEdgeIt.
278 278
    template <typename CutMap>
279 279
    Value minCutMap(const Node& s, ///< The base node.
280 280
                    const Node& t,
281 281
                    ///< The node you want to separate from node \c s.
282 282
                    CutMap& cutMap
283 283
                    ///< The cut will be returned in this map.
284 284
                    /// It must be a \c bool (or convertible) 
285 285
                    /// \ref concepts::ReadWriteMap "ReadWriteMap"
286 286
                    /// on the graph nodes.
287 287
                    ) const {
288 288
      Node sn = s, tn = t;
289 289
      bool s_root=false;
290 290
      Node rn = INVALID;
291 291
      Value value = std::numeric_limits<Value>::max();
292 292
      
293 293
      while (sn != tn) {
294 294
	if ((*_order)[sn] < (*_order)[tn]) {
295 295
	  if ((*_weight)[tn] <= value) {
296 296
	    rn = tn;
297 297
            s_root = false;
298 298
	    value = (*_weight)[tn];
299 299
	  }
300 300
	  tn = (*_pred)[tn];
301 301
	} else {
302 302
	  if ((*_weight)[sn] <= value) {
303 303
	    rn = sn;
304 304
            s_root = true;
305 305
	    value = (*_weight)[sn];
306 306
	  }
307 307
	  sn = (*_pred)[sn];
308 308
	}
309 309
      }
310 310

	
311 311
      typename Graph::template NodeMap<bool> reached(_graph, false);
312
      reached.set(_root, true);
312
      reached[_root] = true;
313 313
      cutMap.set(_root, !s_root);
314
      reached.set(rn, true);
314
      reached[rn] = true;
315 315
      cutMap.set(rn, s_root);
316 316

	
317 317
      std::vector<Node> st;
318 318
      for (NodeIt n(_graph); n != INVALID; ++n) {
319 319
	st.clear();
320 320
        Node nn = n;
321 321
	while (!reached[nn]) {
322 322
	  st.push_back(nn);
323 323
	  nn = (*_pred)[nn];
324 324
	}
325 325
	while (!st.empty()) {
326 326
	  cutMap.set(st.back(), cutMap[nn]);
327 327
	  st.pop_back();
328 328
	}
329 329
      }
330 330
      
331 331
      return value;
332 332
    }
333 333

	
334 334
    ///@}
335 335

	
336 336
    friend class MinCutNodeIt;
337 337

	
338 338
    /// Iterate on the nodes of a minimum cut
339 339
    
340 340
    /// This iterator class lists the nodes of a minimum cut found by
341 341
    /// GomoryHu. Before using it, you must allocate a GomoryHu class,
342 342
    /// and call its \ref GomoryHu::run() "run()" method.
343 343
    ///
344 344
    /// This example counts the nodes in the minimum cut separating \c s from
345 345
    /// \c t.
346 346
    /// \code
347 347
    /// GomoruHu<Graph> gom(g, capacities);
348 348
    /// gom.run();
349 349
    /// int cnt=0;
350 350
    /// for(GomoruHu<Graph>::MinCutNodeIt n(gom,s,t); n!=INVALID; ++n) ++cnt;
351 351
    /// \endcode
352 352
    class MinCutNodeIt
353 353
    {
354 354
      bool _side;
355 355
      typename Graph::NodeIt _node_it;
356 356
      typename Graph::template NodeMap<bool> _cut;
357 357
    public:
358 358
      /// Constructor
359 359

	
360 360
      /// Constructor.
361 361
      ///
362 362
      MinCutNodeIt(GomoryHu const &gomory,
Show white space 96 line context
... ...
@@ -116,732 +116,732 @@
116 116
    typedef typename Digraph::template NodeMap<bool> MinCutMap;
117 117
    MinCutMap* _min_cut_map;
118 118

	
119 119
    Tolerance _tolerance;
120 120

	
121 121
  public:
122 122

	
123 123
    /// \brief Constructor
124 124
    ///
125 125
    /// Constructor of the algorithm class.
126 126
    HaoOrlin(const Digraph& graph, const CapacityMap& capacity,
127 127
             const Tolerance& tolerance = Tolerance()) :
128 128
      _graph(graph), _capacity(&capacity), _flow(0), _source(),
129 129
      _node_num(), _first(), _last(), _next(0), _prev(0),
130 130
      _active(0), _bucket(0), _dormant(), _sets(), _highest(),
131 131
      _excess(0), _source_set(0), _min_cut(), _min_cut_map(0),
132 132
      _tolerance(tolerance) {}
133 133

	
134 134
    ~HaoOrlin() {
135 135
      if (_min_cut_map) {
136 136
        delete _min_cut_map;
137 137
      }
138 138
      if (_source_set) {
139 139
        delete _source_set;
140 140
      }
141 141
      if (_excess) {
142 142
        delete _excess;
143 143
      }
144 144
      if (_next) {
145 145
        delete _next;
146 146
      }
147 147
      if (_prev) {
148 148
        delete _prev;
149 149
      }
150 150
      if (_active) {
151 151
        delete _active;
152 152
      }
153 153
      if (_bucket) {
154 154
        delete _bucket;
155 155
      }
156 156
      if (_flow) {
157 157
        delete _flow;
158 158
      }
159 159
    }
160 160

	
161 161
  private:
162 162

	
163 163
    void activate(const Node& i) {
164
      _active->set(i, true);
164
      (*_active)[i] = true;
165 165

	
166 166
      int bucket = (*_bucket)[i];
167 167

	
168 168
      if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return;
169 169
      //unlace
170
      _next->set((*_prev)[i], (*_next)[i]);
170
      (*_next)[(*_prev)[i]] = (*_next)[i];
171 171
      if ((*_next)[i] != INVALID) {
172
        _prev->set((*_next)[i], (*_prev)[i]);
172
        (*_prev)[(*_next)[i]] = (*_prev)[i];
173 173
      } else {
174 174
        _last[bucket] = (*_prev)[i];
175 175
      }
176 176
      //lace
177
      _next->set(i, _first[bucket]);
178
      _prev->set(_first[bucket], i);
179
      _prev->set(i, INVALID);
177
      (*_next)[i] = _first[bucket];
178
      (*_prev)[_first[bucket]] = i;
179
      (*_prev)[i] = INVALID;
180 180
      _first[bucket] = i;
181 181
    }
182 182

	
183 183
    void deactivate(const Node& i) {
184
      _active->set(i, false);
184
      (*_active)[i] = false;
185 185
      int bucket = (*_bucket)[i];
186 186

	
187 187
      if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return;
188 188

	
189 189
      //unlace
190
      _prev->set((*_next)[i], (*_prev)[i]);
190
      (*_prev)[(*_next)[i]] = (*_prev)[i];
191 191
      if ((*_prev)[i] != INVALID) {
192
        _next->set((*_prev)[i], (*_next)[i]);
192
        (*_next)[(*_prev)[i]] = (*_next)[i];
193 193
      } else {
194 194
        _first[bucket] = (*_next)[i];
195 195
      }
196 196
      //lace
197
      _prev->set(i, _last[bucket]);
198
      _next->set(_last[bucket], i);
199
      _next->set(i, INVALID);
197
      (*_prev)[i] = _last[bucket];
198
      (*_next)[_last[bucket]] = i;
199
      (*_next)[i] = INVALID;
200 200
      _last[bucket] = i;
201 201
    }
202 202

	
203 203
    void addItem(const Node& i, int bucket) {
204 204
      (*_bucket)[i] = bucket;
205 205
      if (_last[bucket] != INVALID) {
206
        _prev->set(i, _last[bucket]);
207
        _next->set(_last[bucket], i);
208
        _next->set(i, INVALID);
206
        (*_prev)[i] = _last[bucket];
207
        (*_next)[_last[bucket]] = i;
208
        (*_next)[i] = INVALID;
209 209
        _last[bucket] = i;
210 210
      } else {
211
        _prev->set(i, INVALID);
211
        (*_prev)[i] = INVALID;
212 212
        _first[bucket] = i;
213
        _next->set(i, INVALID);
213
        (*_next)[i] = INVALID;
214 214
        _last[bucket] = i;
215 215
      }
216 216
    }
217 217

	
218 218
    void findMinCutOut() {
219 219

	
220 220
      for (NodeIt n(_graph); n != INVALID; ++n) {
221
        _excess->set(n, 0);
221
        (*_excess)[n] = 0;
222 222
      }
223 223

	
224 224
      for (ArcIt a(_graph); a != INVALID; ++a) {
225
        _flow->set(a, 0);
225
        (*_flow)[a] = 0;
226 226
      }
227 227

	
228 228
      int bucket_num = 0;
229 229
      std::vector<Node> queue(_node_num);
230 230
      int qfirst = 0, qlast = 0, qsep = 0;
231 231

	
232 232
      {
233 233
        typename Digraph::template NodeMap<bool> reached(_graph, false);
234 234

	
235
        reached.set(_source, true);
235
        reached[_source] = true;
236 236
        bool first_set = true;
237 237

	
238 238
        for (NodeIt t(_graph); t != INVALID; ++t) {
239 239
          if (reached[t]) continue;
240 240
          _sets.push_front(std::list<int>());
241 241

	
242 242
          queue[qlast++] = t;
243
          reached.set(t, true);
243
          reached[t] = true;
244 244

	
245 245
          while (qfirst != qlast) {
246 246
            if (qsep == qfirst) {
247 247
              ++bucket_num;
248 248
              _sets.front().push_front(bucket_num);
249 249
              _dormant[bucket_num] = !first_set;
250 250
              _first[bucket_num] = _last[bucket_num] = INVALID;
251 251
              qsep = qlast;
252 252
            }
253 253

	
254 254
            Node n = queue[qfirst++];
255 255
            addItem(n, bucket_num);
256 256

	
257 257
            for (InArcIt a(_graph, n); a != INVALID; ++a) {
258 258
              Node u = _graph.source(a);
259 259
              if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
260
                reached.set(u, true);
260
                reached[u] = true;
261 261
                queue[qlast++] = u;
262 262
              }
263 263
            }
264 264
          }
265 265
          first_set = false;
266 266
        }
267 267

	
268 268
        ++bucket_num;
269
        _bucket->set(_source, 0);
269
        (*_bucket)[_source] = 0;
270 270
        _dormant[0] = true;
271 271
      }
272
      _source_set->set(_source, true);
272
      (*_source_set)[_source] = true;
273 273

	
274 274
      Node target = _last[_sets.back().back()];
275 275
      {
276 276
        for (OutArcIt a(_graph, _source); a != INVALID; ++a) {
277 277
          if (_tolerance.positive((*_capacity)[a])) {
278 278
            Node u = _graph.target(a);
279
            _flow->set(a, (*_capacity)[a]);
280
            _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
279
            (*_flow)[a] = (*_capacity)[a];
280
            (*_excess)[u] += (*_capacity)[a];
281 281
            if (!(*_active)[u] && u != _source) {
282 282
              activate(u);
283 283
            }
284 284
          }
285 285
        }
286 286

	
287 287
        if ((*_active)[target]) {
288 288
          deactivate(target);
289 289
        }
290 290

	
291 291
        _highest = _sets.back().begin();
292 292
        while (_highest != _sets.back().end() &&
293 293
               !(*_active)[_first[*_highest]]) {
294 294
          ++_highest;
295 295
        }
296 296
      }
297 297

	
298 298
      while (true) {
299 299
        while (_highest != _sets.back().end()) {
300 300
          Node n = _first[*_highest];
301 301
          Value excess = (*_excess)[n];
302 302
          int next_bucket = _node_num;
303 303

	
304 304
          int under_bucket;
305 305
          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
306 306
            under_bucket = -1;
307 307
          } else {
308 308
            under_bucket = *(++std::list<int>::iterator(_highest));
309 309
          }
310 310

	
311 311
          for (OutArcIt a(_graph, n); a != INVALID; ++a) {
312 312
            Node v = _graph.target(a);
313 313
            if (_dormant[(*_bucket)[v]]) continue;
314 314
            Value rem = (*_capacity)[a] - (*_flow)[a];
315 315
            if (!_tolerance.positive(rem)) continue;
316 316
            if ((*_bucket)[v] == under_bucket) {
317 317
              if (!(*_active)[v] && v != target) {
318 318
                activate(v);
319 319
              }
320 320
              if (!_tolerance.less(rem, excess)) {
321
                _flow->set(a, (*_flow)[a] + excess);
322
                _excess->set(v, (*_excess)[v] + excess);
321
                (*_flow)[a] += excess;
322
                (*_excess)[v] += excess;
323 323
                excess = 0;
324 324
                goto no_more_push;
325 325
              } else {
326 326
                excess -= rem;
327
                _excess->set(v, (*_excess)[v] + rem);
328
                _flow->set(a, (*_capacity)[a]);
327
                (*_excess)[v] += rem;
328
                (*_flow)[a] = (*_capacity)[a];
329 329
              }
330 330
            } else if (next_bucket > (*_bucket)[v]) {
331 331
              next_bucket = (*_bucket)[v];
332 332
            }
333 333
          }
334 334

	
335 335
          for (InArcIt a(_graph, n); a != INVALID; ++a) {
336 336
            Node v = _graph.source(a);
337 337
            if (_dormant[(*_bucket)[v]]) continue;
338 338
            Value rem = (*_flow)[a];
339 339
            if (!_tolerance.positive(rem)) continue;
340 340
            if ((*_bucket)[v] == under_bucket) {
341 341
              if (!(*_active)[v] && v != target) {
342 342
                activate(v);
343 343
              }
344 344
              if (!_tolerance.less(rem, excess)) {
345
                _flow->set(a, (*_flow)[a] - excess);
346
                _excess->set(v, (*_excess)[v] + excess);
345
                (*_flow)[a] -= excess;
346
                (*_excess)[v] += excess;
347 347
                excess = 0;
348 348
                goto no_more_push;
349 349
              } else {
350 350
                excess -= rem;
351
                _excess->set(v, (*_excess)[v] + rem);
352
                _flow->set(a, 0);
351
                (*_excess)[v] += rem;
352
                (*_flow)[a] = 0;
353 353
              }
354 354
            } else if (next_bucket > (*_bucket)[v]) {
355 355
              next_bucket = (*_bucket)[v];
356 356
            }
357 357
          }
358 358

	
359 359
        no_more_push:
360 360

	
361
          _excess->set(n, excess);
361
          (*_excess)[n] = excess;
362 362

	
363 363
          if (excess != 0) {
364 364
            if ((*_next)[n] == INVALID) {
365 365
              typename std::list<std::list<int> >::iterator new_set =
366 366
                _sets.insert(--_sets.end(), std::list<int>());
367 367
              new_set->splice(new_set->end(), _sets.back(),
368 368
                              _sets.back().begin(), ++_highest);
369 369
              for (std::list<int>::iterator it = new_set->begin();
370 370
                   it != new_set->end(); ++it) {
371 371
                _dormant[*it] = true;
372 372
              }
373 373
              while (_highest != _sets.back().end() &&
374 374
                     !(*_active)[_first[*_highest]]) {
375 375
                ++_highest;
376 376
              }
377 377
            } else if (next_bucket == _node_num) {
378 378
              _first[(*_bucket)[n]] = (*_next)[n];
379
              _prev->set((*_next)[n], INVALID);
379
              (*_prev)[(*_next)[n]] = INVALID;
380 380

	
381 381
              std::list<std::list<int> >::iterator new_set =
382 382
                _sets.insert(--_sets.end(), std::list<int>());
383 383

	
384 384
              new_set->push_front(bucket_num);
385
              _bucket->set(n, bucket_num);
385
              (*_bucket)[n] = bucket_num;
386 386
              _first[bucket_num] = _last[bucket_num] = n;
387
              _next->set(n, INVALID);
388
              _prev->set(n, INVALID);
387
              (*_next)[n] = INVALID;
388
              (*_prev)[n] = INVALID;
389 389
              _dormant[bucket_num] = true;
390 390
              ++bucket_num;
391 391

	
392 392
              while (_highest != _sets.back().end() &&
393 393
                     !(*_active)[_first[*_highest]]) {
394 394
                ++_highest;
395 395
              }
396 396
            } else {
397 397
              _first[*_highest] = (*_next)[n];
398
              _prev->set((*_next)[n], INVALID);
398
              (*_prev)[(*_next)[n]] = INVALID;
399 399

	
400 400
              while (next_bucket != *_highest) {
401 401
                --_highest;
402 402
              }
403 403

	
404 404
              if (_highest == _sets.back().begin()) {
405 405
                _sets.back().push_front(bucket_num);
406 406
                _dormant[bucket_num] = false;
407 407
                _first[bucket_num] = _last[bucket_num] = INVALID;
408 408
                ++bucket_num;
409 409
              }
410 410
              --_highest;
411 411

	
412
              _bucket->set(n, *_highest);
413
              _next->set(n, _first[*_highest]);
412
              (*_bucket)[n] = *_highest;
413
              (*_next)[n] = _first[*_highest];
414 414
              if (_first[*_highest] != INVALID) {
415
                _prev->set(_first[*_highest], n);
415
                (*_prev)[_first[*_highest]] = n;
416 416
              } else {
417 417
                _last[*_highest] = n;
418 418
              }
419 419
              _first[*_highest] = n;
420 420
            }
421 421
          } else {
422 422

	
423 423
            deactivate(n);
424 424
            if (!(*_active)[_first[*_highest]]) {
425 425
              ++_highest;
426 426
              if (_highest != _sets.back().end() &&
427 427
                  !(*_active)[_first[*_highest]]) {
428 428
                _highest = _sets.back().end();
429 429
              }
430 430
            }
431 431
          }
432 432
        }
433 433

	
434 434
        if ((*_excess)[target] < _min_cut) {
435 435
          _min_cut = (*_excess)[target];
436 436
          for (NodeIt i(_graph); i != INVALID; ++i) {
437
            _min_cut_map->set(i, true);
437
            (*_min_cut_map)[i] = true;
438 438
          }
439 439
          for (std::list<int>::iterator it = _sets.back().begin();
440 440
               it != _sets.back().end(); ++it) {
441 441
            Node n = _first[*it];
442 442
            while (n != INVALID) {
443
              _min_cut_map->set(n, false);
443
              (*_min_cut_map)[n] = false;
444 444
              n = (*_next)[n];
445 445
            }
446 446
          }
447 447
        }
448 448

	
449 449
        {
450 450
          Node new_target;
451 451
          if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
452 452
            if ((*_next)[target] == INVALID) {
453 453
              _last[(*_bucket)[target]] = (*_prev)[target];
454 454
              new_target = (*_prev)[target];
455 455
            } else {
456
              _prev->set((*_next)[target], (*_prev)[target]);
456
              (*_prev)[(*_next)[target]] = (*_prev)[target];
457 457
              new_target = (*_next)[target];
458 458
            }
459 459
            if ((*_prev)[target] == INVALID) {
460 460
              _first[(*_bucket)[target]] = (*_next)[target];
461 461
            } else {
462
              _next->set((*_prev)[target], (*_next)[target]);
462
              (*_next)[(*_prev)[target]] = (*_next)[target];
463 463
            }
464 464
          } else {
465 465
            _sets.back().pop_back();
466 466
            if (_sets.back().empty()) {
467 467
              _sets.pop_back();
468 468
              if (_sets.empty())
469 469
                break;
470 470
              for (std::list<int>::iterator it = _sets.back().begin();
471 471
                   it != _sets.back().end(); ++it) {
472 472
                _dormant[*it] = false;
473 473
              }
474 474
            }
475 475
            new_target = _last[_sets.back().back()];
476 476
          }
477 477

	
478
          _bucket->set(target, 0);
478
          (*_bucket)[target] = 0;
479 479

	
480
          _source_set->set(target, true);
480
          (*_source_set)[target] = true;
481 481
          for (OutArcIt a(_graph, target); a != INVALID; ++a) {
482 482
            Value rem = (*_capacity)[a] - (*_flow)[a];
483 483
            if (!_tolerance.positive(rem)) continue;
484 484
            Node v = _graph.target(a);
485 485
            if (!(*_active)[v] && !(*_source_set)[v]) {
486 486
              activate(v);
487 487
            }
488
            _excess->set(v, (*_excess)[v] + rem);
489
            _flow->set(a, (*_capacity)[a]);
488
            (*_excess)[v] += rem;
489
            (*_flow)[a] = (*_capacity)[a];
490 490
          }
491 491

	
492 492
          for (InArcIt a(_graph, target); a != INVALID; ++a) {
493 493
            Value rem = (*_flow)[a];
494 494
            if (!_tolerance.positive(rem)) continue;
495 495
            Node v = _graph.source(a);
496 496
            if (!(*_active)[v] && !(*_source_set)[v]) {
497 497
              activate(v);
498 498
            }
499
            _excess->set(v, (*_excess)[v] + rem);
500
            _flow->set(a, 0);
499
            (*_excess)[v] += rem;
500
            (*_flow)[a] = 0;
501 501
          }
502 502

	
503 503
          target = new_target;
504 504
          if ((*_active)[target]) {
505 505
            deactivate(target);
506 506
          }
507 507

	
508 508
          _highest = _sets.back().begin();
509 509
          while (_highest != _sets.back().end() &&
510 510
                 !(*_active)[_first[*_highest]]) {
511 511
            ++_highest;
512 512
          }
513 513
        }
514 514
      }
515 515
    }
516 516

	
517 517
    void findMinCutIn() {
518 518

	
519 519
      for (NodeIt n(_graph); n != INVALID; ++n) {
520
        _excess->set(n, 0);
520
        (*_excess)[n] = 0;
521 521
      }
522 522

	
523 523
      for (ArcIt a(_graph); a != INVALID; ++a) {
524
        _flow->set(a, 0);
524
        (*_flow)[a] = 0;
525 525
      }
526 526

	
527 527
      int bucket_num = 0;
528 528
      std::vector<Node> queue(_node_num);
529 529
      int qfirst = 0, qlast = 0, qsep = 0;
530 530

	
531 531
      {
532 532
        typename Digraph::template NodeMap<bool> reached(_graph, false);
533 533

	
534
        reached.set(_source, true);
534
        reached[_source] = true;
535 535

	
536 536
        bool first_set = true;
537 537

	
538 538
        for (NodeIt t(_graph); t != INVALID; ++t) {
539 539
          if (reached[t]) continue;
540 540
          _sets.push_front(std::list<int>());
541 541

	
542 542
          queue[qlast++] = t;
543
          reached.set(t, true);
543
          reached[t] = true;
544 544

	
545 545
          while (qfirst != qlast) {
546 546
            if (qsep == qfirst) {
547 547
              ++bucket_num;
548 548
              _sets.front().push_front(bucket_num);
549 549
              _dormant[bucket_num] = !first_set;
550 550
              _first[bucket_num] = _last[bucket_num] = INVALID;
551 551
              qsep = qlast;
552 552
            }
553 553

	
554 554
            Node n = queue[qfirst++];
555 555
            addItem(n, bucket_num);
556 556

	
557 557
            for (OutArcIt a(_graph, n); a != INVALID; ++a) {
558 558
              Node u = _graph.target(a);
559 559
              if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
560
                reached.set(u, true);
560
                reached[u] = true;
561 561
                queue[qlast++] = u;
562 562
              }
563 563
            }
564 564
          }
565 565
          first_set = false;
566 566
        }
567 567

	
568 568
        ++bucket_num;
569
        _bucket->set(_source, 0);
569
        (*_bucket)[_source] = 0;
570 570
        _dormant[0] = true;
571 571
      }
572
      _source_set->set(_source, true);
572
      (*_source_set)[_source] = true;
573 573

	
574 574
      Node target = _last[_sets.back().back()];
575 575
      {
576 576
        for (InArcIt a(_graph, _source); a != INVALID; ++a) {
577 577
          if (_tolerance.positive((*_capacity)[a])) {
578 578
            Node u = _graph.source(a);
579
            _flow->set(a, (*_capacity)[a]);
580
            _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
579
            (*_flow)[a] = (*_capacity)[a];
580
            (*_excess)[u] += (*_capacity)[a];
581 581
            if (!(*_active)[u] && u != _source) {
582 582
              activate(u);
583 583
            }
584 584
          }
585 585
        }
586 586
        if ((*_active)[target]) {
587 587
          deactivate(target);
588 588
        }
589 589

	
590 590
        _highest = _sets.back().begin();
591 591
        while (_highest != _sets.back().end() &&
592 592
               !(*_active)[_first[*_highest]]) {
593 593
          ++_highest;
594 594
        }
595 595
      }
596 596

	
597 597

	
598 598
      while (true) {
599 599
        while (_highest != _sets.back().end()) {
600 600
          Node n = _first[*_highest];
601 601
          Value excess = (*_excess)[n];
602 602
          int next_bucket = _node_num;
603 603

	
604 604
          int under_bucket;
605 605
          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
606 606
            under_bucket = -1;
607 607
          } else {
608 608
            under_bucket = *(++std::list<int>::iterator(_highest));
609 609
          }
610 610

	
611 611
          for (InArcIt a(_graph, n); a != INVALID; ++a) {
612 612
            Node v = _graph.source(a);
613 613
            if (_dormant[(*_bucket)[v]]) continue;
614 614
            Value rem = (*_capacity)[a] - (*_flow)[a];
615 615
            if (!_tolerance.positive(rem)) continue;
616 616
            if ((*_bucket)[v] == under_bucket) {
617 617
              if (!(*_active)[v] && v != target) {
618 618
                activate(v);
619 619
              }
620 620
              if (!_tolerance.less(rem, excess)) {
621
                _flow->set(a, (*_flow)[a] + excess);
622
                _excess->set(v, (*_excess)[v] + excess);
621
                (*_flow)[a] += excess;
622
                (*_excess)[v] += excess;
623 623
                excess = 0;
624 624
                goto no_more_push;
625 625
              } else {
626 626
                excess -= rem;
627
                _excess->set(v, (*_excess)[v] + rem);
628
                _flow->set(a, (*_capacity)[a]);
627
                (*_excess)[v] += rem;
628
                (*_flow)[a] = (*_capacity)[a];
629 629
              }
630 630
            } else if (next_bucket > (*_bucket)[v]) {
631 631
              next_bucket = (*_bucket)[v];
632 632
            }
633 633
          }
634 634

	
635 635
          for (OutArcIt a(_graph, n); a != INVALID; ++a) {
636 636
            Node v = _graph.target(a);
637 637
            if (_dormant[(*_bucket)[v]]) continue;
638 638
            Value rem = (*_flow)[a];
639 639
            if (!_tolerance.positive(rem)) continue;
640 640
            if ((*_bucket)[v] == under_bucket) {
641 641
              if (!(*_active)[v] && v != target) {
642 642
                activate(v);
643 643
              }
644 644
              if (!_tolerance.less(rem, excess)) {
645
                _flow->set(a, (*_flow)[a] - excess);
646
                _excess->set(v, (*_excess)[v] + excess);
645
                (*_flow)[a] -= excess;
646
                (*_excess)[v] += excess;
647 647
                excess = 0;
648 648
                goto no_more_push;
649 649
              } else {
650 650
                excess -= rem;
651
                _excess->set(v, (*_excess)[v] + rem);
652
                _flow->set(a, 0);
651
                (*_excess)[v] += rem;
652
                (*_flow)[a] = 0;
653 653
              }
654 654
            } else if (next_bucket > (*_bucket)[v]) {
655 655
              next_bucket = (*_bucket)[v];
656 656
            }
657 657
          }
658 658

	
659 659
        no_more_push:
660 660

	
661
          _excess->set(n, excess);
661
          (*_excess)[n] = excess;
662 662

	
663 663
          if (excess != 0) {
664 664
            if ((*_next)[n] == INVALID) {
665 665
              typename std::list<std::list<int> >::iterator new_set =
666 666
                _sets.insert(--_sets.end(), std::list<int>());
667 667
              new_set->splice(new_set->end(), _sets.back(),
668 668
                              _sets.back().begin(), ++_highest);
669 669
              for (std::list<int>::iterator it = new_set->begin();
670 670
                   it != new_set->end(); ++it) {
671 671
                _dormant[*it] = true;
672 672
              }
673 673
              while (_highest != _sets.back().end() &&
674 674
                     !(*_active)[_first[*_highest]]) {
675 675
                ++_highest;
676 676
              }
677 677
            } else if (next_bucket == _node_num) {
678 678
              _first[(*_bucket)[n]] = (*_next)[n];
679
              _prev->set((*_next)[n], INVALID);
679
              (*_prev)[(*_next)[n]] = INVALID;
680 680

	
681 681
              std::list<std::list<int> >::iterator new_set =
682 682
                _sets.insert(--_sets.end(), std::list<int>());
683 683

	
684 684
              new_set->push_front(bucket_num);
685
              _bucket->set(n, bucket_num);
685
              (*_bucket)[n] = bucket_num;
686 686
              _first[bucket_num] = _last[bucket_num] = n;
687
              _next->set(n, INVALID);
688
              _prev->set(n, INVALID);
687
              (*_next)[n] = INVALID;
688
              (*_prev)[n] = INVALID;
689 689
              _dormant[bucket_num] = true;
690 690
              ++bucket_num;
691 691

	
692 692
              while (_highest != _sets.back().end() &&
693 693
                     !(*_active)[_first[*_highest]]) {
694 694
                ++_highest;
695 695
              }
696 696
            } else {
697 697
              _first[*_highest] = (*_next)[n];
698
              _prev->set((*_next)[n], INVALID);
698
              (*_prev)[(*_next)[n]] = INVALID;
699 699

	
700 700
              while (next_bucket != *_highest) {
701 701
                --_highest;
702 702
              }
703 703
              if (_highest == _sets.back().begin()) {
704 704
                _sets.back().push_front(bucket_num);
705 705
                _dormant[bucket_num] = false;
706 706
                _first[bucket_num] = _last[bucket_num] = INVALID;
707 707
                ++bucket_num;
708 708
              }
709 709
              --_highest;
710 710

	
711
              _bucket->set(n, *_highest);
712
              _next->set(n, _first[*_highest]);
711
              (*_bucket)[n] = *_highest;
712
              (*_next)[n] = _first[*_highest];
713 713
              if (_first[*_highest] != INVALID) {
714
                _prev->set(_first[*_highest], n);
714
                (*_prev)[_first[*_highest]] = n;
715 715
              } else {
716 716
                _last[*_highest] = n;
717 717
              }
718 718
              _first[*_highest] = n;
719 719
            }
720 720
          } else {
721 721

	
722 722
            deactivate(n);
723 723
            if (!(*_active)[_first[*_highest]]) {
724 724
              ++_highest;
725 725
              if (_highest != _sets.back().end() &&
726 726
                  !(*_active)[_first[*_highest]]) {
727 727
                _highest = _sets.back().end();
728 728
              }
729 729
            }
730 730
          }
731 731
        }
732 732

	
733 733
        if ((*_excess)[target] < _min_cut) {
734 734
          _min_cut = (*_excess)[target];
735 735
          for (NodeIt i(_graph); i != INVALID; ++i) {
736
            _min_cut_map->set(i, false);
736
            (*_min_cut_map)[i] = false;
737 737
          }
738 738
          for (std::list<int>::iterator it = _sets.back().begin();
739 739
               it != _sets.back().end(); ++it) {
740 740
            Node n = _first[*it];
741 741
            while (n != INVALID) {
742
              _min_cut_map->set(n, true);
742
              (*_min_cut_map)[n] = true;
743 743
              n = (*_next)[n];
744 744
            }
745 745
          }
746 746
        }
747 747

	
748 748
        {
749 749
          Node new_target;
750 750
          if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
751 751
            if ((*_next)[target] == INVALID) {
752 752
              _last[(*_bucket)[target]] = (*_prev)[target];
753 753
              new_target = (*_prev)[target];
754 754
            } else {
755
              _prev->set((*_next)[target], (*_prev)[target]);
755
              (*_prev)[(*_next)[target]] = (*_prev)[target];
756 756
              new_target = (*_next)[target];
757 757
            }
758 758
            if ((*_prev)[target] == INVALID) {
759 759
              _first[(*_bucket)[target]] = (*_next)[target];
760 760
            } else {
761
              _next->set((*_prev)[target], (*_next)[target]);
761
              (*_next)[(*_prev)[target]] = (*_next)[target];
762 762
            }
763 763
          } else {
764 764
            _sets.back().pop_back();
765 765
            if (_sets.back().empty()) {
766 766
              _sets.pop_back();
767 767
              if (_sets.empty())
768 768
                break;
769 769
              for (std::list<int>::iterator it = _sets.back().begin();
770 770
                   it != _sets.back().end(); ++it) {
771 771
                _dormant[*it] = false;
772 772
              }
773 773
            }
774 774
            new_target = _last[_sets.back().back()];
775 775
          }
776 776

	
777
          _bucket->set(target, 0);
777
          (*_bucket)[target] = 0;
778 778

	
779
          _source_set->set(target, true);
779
          (*_source_set)[target] = true;
780 780
          for (InArcIt a(_graph, target); a != INVALID; ++a) {
781 781
            Value rem = (*_capacity)[a] - (*_flow)[a];
782 782
            if (!_tolerance.positive(rem)) continue;
783 783
            Node v = _graph.source(a);
784 784
            if (!(*_active)[v] && !(*_source_set)[v]) {
785 785
              activate(v);
786 786
            }
787
            _excess->set(v, (*_excess)[v] + rem);
788
            _flow->set(a, (*_capacity)[a]);
787
            (*_excess)[v] += rem;
788
            (*_flow)[a] = (*_capacity)[a];
789 789
          }
790 790

	
791 791
          for (OutArcIt a(_graph, target); a != INVALID; ++a) {
792 792
            Value rem = (*_flow)[a];
793 793
            if (!_tolerance.positive(rem)) continue;
794 794
            Node v = _graph.target(a);
795 795
            if (!(*_active)[v] && !(*_source_set)[v]) {
796 796
              activate(v);
797 797
            }
798
            _excess->set(v, (*_excess)[v] + rem);
799
            _flow->set(a, 0);
798
            (*_excess)[v] += rem;
799
            (*_flow)[a] = 0;
800 800
          }
801 801

	
802 802
          target = new_target;
803 803
          if ((*_active)[target]) {
804 804
            deactivate(target);
805 805
          }
806 806

	
807 807
          _highest = _sets.back().begin();
808 808
          while (_highest != _sets.back().end() &&
809 809
                 !(*_active)[_first[*_highest]]) {
810 810
            ++_highest;
811 811
          }
812 812
        }
813 813
      }
814 814
    }
815 815

	
816 816
  public:
817 817

	
818 818
    /// \name Execution control
819 819
    /// The simplest way to execute the algorithm is to use
820 820
    /// one of the member functions called \ref run().
821 821
    /// \n
822 822
    /// If you need more control on the execution,
823 823
    /// first you must call \ref init(), then the \ref calculateIn() or
824 824
    /// \ref calculateOut() functions.
825 825

	
826 826
    /// @{
827 827

	
828 828
    /// \brief Initializes the internal data structures.
829 829
    ///
830 830
    /// Initializes the internal data structures. It creates
831 831
    /// the maps, residual graph adaptors and some bucket structures
832 832
    /// for the algorithm.
833 833
    void init() {
834 834
      init(NodeIt(_graph));
835 835
    }
836 836

	
837 837
    /// \brief Initializes the internal data structures.
838 838
    ///
839 839
    /// Initializes the internal data structures. It creates
840 840
    /// the maps, residual graph adaptor and some bucket structures
841 841
    /// for the algorithm. Node \c source  is used as the push-relabel
842 842
    /// algorithm's source.
843 843
    void init(const Node& source) {
844 844
      _source = source;
845 845

	
846 846
      _node_num = countNodes(_graph);
847 847

	
Show white space 96 line context
... ...
@@ -237,327 +237,327 @@
237 237
        while (true) {
238 238
          if ((*_matching)[left] == INVALID) break;
239 239
          left = _graph.target((*_matching)[left]);
240 240
          left = (*_blossom_rep)[_blossom_set->
241 241
                                 find(_graph.target((*_ear)[left]))];
242 242
          if (right_set.find(left) != right_set.end()) {
243 243
            nca = left;
244 244
            break;
245 245
          }
246 246
          left_set.insert(left);
247 247

	
248 248
          if ((*_matching)[right] == INVALID) break;
249 249
          right = _graph.target((*_matching)[right]);
250 250
          right = (*_blossom_rep)[_blossom_set->
251 251
                                  find(_graph.target((*_ear)[right]))];
252 252
          if (left_set.find(right) != left_set.end()) {
253 253
            nca = right;
254 254
            break;
255 255
          }
256 256
          right_set.insert(right);
257 257
        }
258 258

	
259 259
        if (nca == INVALID) {
260 260
          if ((*_matching)[left] == INVALID) {
261 261
            nca = right;
262 262
            while (left_set.find(nca) == left_set.end()) {
263 263
              nca = _graph.target((*_matching)[nca]);
264 264
              nca =(*_blossom_rep)[_blossom_set->
265 265
                                   find(_graph.target((*_ear)[nca]))];
266 266
            }
267 267
          } else {
268 268
            nca = left;
269 269
            while (right_set.find(nca) == right_set.end()) {
270 270
              nca = _graph.target((*_matching)[nca]);
271 271
              nca = (*_blossom_rep)[_blossom_set->
272 272
                                   find(_graph.target((*_ear)[nca]))];
273 273
            }
274 274
          }
275 275
        }
276 276
      }
277 277

	
278 278
      {
279 279

	
280 280
        Node node = _graph.u(e);
281 281
        Arc arc = _graph.direct(e, true);
282 282
        Node base = (*_blossom_rep)[_blossom_set->find(node)];
283 283

	
284 284
        while (base != nca) {
285
          _ear->set(node, arc);
285
          (*_ear)[node] = arc;
286 286

	
287 287
          Node n = node;
288 288
          while (n != base) {
289 289
            n = _graph.target((*_matching)[n]);
290 290
            Arc a = (*_ear)[n];
291 291
            n = _graph.target(a);
292
            _ear->set(n, _graph.oppositeArc(a));
292
            (*_ear)[n] = _graph.oppositeArc(a);
293 293
          }
294 294
          node = _graph.target((*_matching)[base]);
295 295
          _tree_set->erase(base);
296 296
          _tree_set->erase(node);
297 297
          _blossom_set->insert(node, _blossom_set->find(base));
298
          _status->set(node, EVEN);
298
          (*_status)[node] = EVEN;
299 299
          _node_queue[_last++] = node;
300 300
          arc = _graph.oppositeArc((*_ear)[node]);
301 301
          node = _graph.target((*_ear)[node]);
302 302
          base = (*_blossom_rep)[_blossom_set->find(node)];
303 303
          _blossom_set->join(_graph.target(arc), base);
304 304
        }
305 305
      }
306 306

	
307
      _blossom_rep->set(_blossom_set->find(nca), nca);
307
      (*_blossom_rep)[_blossom_set->find(nca)] = nca;
308 308

	
309 309
      {
310 310

	
311 311
        Node node = _graph.v(e);
312 312
        Arc arc = _graph.direct(e, false);
313 313
        Node base = (*_blossom_rep)[_blossom_set->find(node)];
314 314

	
315 315
        while (base != nca) {
316
          _ear->set(node, arc);
316
          (*_ear)[node] = arc;
317 317

	
318 318
          Node n = node;
319 319
          while (n != base) {
320 320
            n = _graph.target((*_matching)[n]);
321 321
            Arc a = (*_ear)[n];
322 322
            n = _graph.target(a);
323
            _ear->set(n, _graph.oppositeArc(a));
323
            (*_ear)[n] = _graph.oppositeArc(a);
324 324
          }
325 325
          node = _graph.target((*_matching)[base]);
326 326
          _tree_set->erase(base);
327 327
          _tree_set->erase(node);
328 328
          _blossom_set->insert(node, _blossom_set->find(base));
329
          _status->set(node, EVEN);
329
          (*_status)[node] = EVEN;
330 330
          _node_queue[_last++] = node;
331 331
          arc = _graph.oppositeArc((*_ear)[node]);
332 332
          node = _graph.target((*_ear)[node]);
333 333
          base = (*_blossom_rep)[_blossom_set->find(node)];
334 334
          _blossom_set->join(_graph.target(arc), base);
335 335
        }
336 336
      }
337 337

	
338
      _blossom_rep->set(_blossom_set->find(nca), nca);
338
      (*_blossom_rep)[_blossom_set->find(nca)] = nca;
339 339
    }
340 340

	
341 341

	
342 342

	
343 343
    void extendOnArc(const Arc& a) {
344 344
      Node base = _graph.source(a);
345 345
      Node odd = _graph.target(a);
346 346

	
347
      _ear->set(odd, _graph.oppositeArc(a));
347
      (*_ear)[odd] = _graph.oppositeArc(a);
348 348
      Node even = _graph.target((*_matching)[odd]);
349
      _blossom_rep->set(_blossom_set->insert(even), even);
350
      _status->set(odd, ODD);
351
      _status->set(even, EVEN);
349
      (*_blossom_rep)[_blossom_set->insert(even)] = even;
350
      (*_status)[odd] = ODD;
351
      (*_status)[even] = EVEN;
352 352
      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
353 353
      _tree_set->insert(odd, tree);
354 354
      _tree_set->insert(even, tree);
355 355
      _node_queue[_last++] = even;
356 356

	
357 357
    }
358 358

	
359 359
    void augmentOnArc(const Arc& a) {
360 360
      Node even = _graph.source(a);
361 361
      Node odd = _graph.target(a);
362 362

	
363 363
      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
364 364

	
365
      _matching->set(odd, _graph.oppositeArc(a));
366
      _status->set(odd, MATCHED);
365
      (*_matching)[odd] = _graph.oppositeArc(a);
366
      (*_status)[odd] = MATCHED;
367 367

	
368 368
      Arc arc = (*_matching)[even];
369
      _matching->set(even, a);
369
      (*_matching)[even] = a;
370 370

	
371 371
      while (arc != INVALID) {
372 372
        odd = _graph.target(arc);
373 373
        arc = (*_ear)[odd];
374 374
        even = _graph.target(arc);
375
        _matching->set(odd, arc);
375
        (*_matching)[odd] = arc;
376 376
        arc = (*_matching)[even];
377
        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
377
        (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]);
378 378
      }
379 379

	
380 380
      for (typename TreeSet::ItemIt it(*_tree_set, tree);
381 381
           it != INVALID; ++it) {
382 382
        if ((*_status)[it] == ODD) {
383
          _status->set(it, MATCHED);
383
          (*_status)[it] = MATCHED;
384 384
        } else {
385 385
          int blossom = _blossom_set->find(it);
386 386
          for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
387 387
               jt != INVALID; ++jt) {
388
            _status->set(jt, MATCHED);
388
            (*_status)[jt] = MATCHED;
389 389
          }
390 390
          _blossom_set->eraseClass(blossom);
391 391
        }
392 392
      }
393 393
      _tree_set->eraseClass(tree);
394 394

	
395 395
    }
396 396

	
397 397
  public:
398 398

	
399 399
    /// \brief Constructor
400 400
    ///
401 401
    /// Constructor.
402 402
    MaxMatching(const Graph& graph)
403 403
      : _graph(graph), _matching(0), _status(0), _ear(0),
404 404
        _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
405 405
        _tree_set_index(0), _tree_set(0) {}
406 406

	
407 407
    ~MaxMatching() {
408 408
      destroyStructures();
409 409
    }
410 410

	
411 411
    /// \name Execution control
412 412
    /// The simplest way to execute the algorithm is to use the
413 413
    /// \c run() member function.
414 414
    /// \n
415 415

	
416 416
    /// If you need better control on the execution, you must call
417 417
    /// \ref init(), \ref greedyInit() or \ref matchingInit()
418 418
    /// functions first, then you can start the algorithm with the \ref
419 419
    /// startSparse() or startDense() functions.
420 420

	
421 421
    ///@{
422 422

	
423 423
    /// \brief Sets the actual matching to the empty matching.
424 424
    ///
425 425
    /// Sets the actual matching to the empty matching.
426 426
    ///
427 427
    void init() {
428 428
      createStructures();
429 429
      for(NodeIt n(_graph); n != INVALID; ++n) {
430
        _matching->set(n, INVALID);
431
        _status->set(n, UNMATCHED);
430
        (*_matching)[n] = INVALID;
431
        (*_status)[n] = UNMATCHED;
432 432
      }
433 433
    }
434 434

	
435 435
    ///\brief Finds an initial matching in a greedy way
436 436
    ///
437 437
    ///It finds an initial matching in a greedy way.
438 438
    void greedyInit() {
439 439
      createStructures();
440 440
      for (NodeIt n(_graph); n != INVALID; ++n) {
441
        _matching->set(n, INVALID);
442
        _status->set(n, UNMATCHED);
441
        (*_matching)[n] = INVALID;
442
        (*_status)[n] = UNMATCHED;
443 443
      }
444 444
      for (NodeIt n(_graph); n != INVALID; ++n) {
445 445
        if ((*_matching)[n] == INVALID) {
446 446
          for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
447 447
            Node v = _graph.target(a);
448 448
            if ((*_matching)[v] == INVALID && v != n) {
449
              _matching->set(n, a);
450
              _status->set(n, MATCHED);
451
              _matching->set(v, _graph.oppositeArc(a));
452
              _status->set(v, MATCHED);
449
              (*_matching)[n] = a;
450
              (*_status)[n] = MATCHED;
451
              (*_matching)[v] = _graph.oppositeArc(a);
452
              (*_status)[v] = MATCHED;
453 453
              break;
454 454
            }
455 455
          }
456 456
        }
457 457
      }
458 458
    }
459 459

	
460 460

	
461 461
    /// \brief Initialize the matching from a map containing.
462 462
    ///
463 463
    /// Initialize the matching from a \c bool valued \c Edge map. This
464 464
    /// map must have the property that there are no two incident edges
465 465
    /// with true value, ie. it contains a matching.
466 466
    /// \return \c true if the map contains a matching.
467 467
    template <typename MatchingMap>
468 468
    bool matchingInit(const MatchingMap& matching) {
469 469
      createStructures();
470 470

	
471 471
      for (NodeIt n(_graph); n != INVALID; ++n) {
472
        _matching->set(n, INVALID);
473
        _status->set(n, UNMATCHED);
472
        (*_matching)[n] = INVALID;
473
        (*_status)[n] = UNMATCHED;
474 474
      }
475 475
      for(EdgeIt e(_graph); e!=INVALID; ++e) {
476 476
        if (matching[e]) {
477 477

	
478 478
          Node u = _graph.u(e);
479 479
          if ((*_matching)[u] != INVALID) return false;
480
          _matching->set(u, _graph.direct(e, true));
481
          _status->set(u, MATCHED);
480
          (*_matching)[u] = _graph.direct(e, true);
481
          (*_status)[u] = MATCHED;
482 482

	
483 483
          Node v = _graph.v(e);
484 484
          if ((*_matching)[v] != INVALID) return false;
485
          _matching->set(v, _graph.direct(e, false));
486
          _status->set(v, MATCHED);
485
          (*_matching)[v] = _graph.direct(e, false);
486
          (*_status)[v] = MATCHED;
487 487
        }
488 488
      }
489 489
      return true;
490 490
    }
491 491

	
492 492
    /// \brief Starts Edmonds' algorithm
493 493
    ///
494 494
    /// If runs the original Edmonds' algorithm.
495 495
    void startSparse() {
496 496
      for(NodeIt n(_graph); n != INVALID; ++n) {
497 497
        if ((*_status)[n] == UNMATCHED) {
498 498
          (*_blossom_rep)[_blossom_set->insert(n)] = n;
499 499
          _tree_set->insert(n);
500
          _status->set(n, EVEN);
500
          (*_status)[n] = EVEN;
501 501
          processSparse(n);
502 502
        }
503 503
      }
504 504
    }
505 505

	
506 506
    /// \brief Starts Edmonds' algorithm.
507 507
    ///
508 508
    /// It runs Edmonds' algorithm with a heuristic of postponing
509 509
    /// shrinks, therefore resulting in a faster algorithm for dense graphs.
510 510
    void startDense() {
511 511
      for(NodeIt n(_graph); n != INVALID; ++n) {
512 512
        if ((*_status)[n] == UNMATCHED) {
513 513
          (*_blossom_rep)[_blossom_set->insert(n)] = n;
514 514
          _tree_set->insert(n);
515
          _status->set(n, EVEN);
515
          (*_status)[n] = EVEN;
516 516
          processDense(n);
517 517
        }
518 518
      }
519 519
    }
520 520

	
521 521

	
522 522
    /// \brief Runs Edmonds' algorithm
523 523
    ///
524 524
    /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
525 525
    /// or Edmonds' algorithm with a heuristic of
526 526
    /// postponing shrinks for dense graphs.
527 527
    void run() {
528 528
      if (countEdges(_graph) < 2 * countNodes(_graph)) {
529 529
        greedyInit();
530 530
        startSparse();
531 531
      } else {
532 532
        init();
533 533
        startDense();
534 534
      }
535 535
    }
536 536

	
537 537
    /// @}
538 538

	
539 539
    /// \name Primal solution
540 540
    /// Functions to get the primal solution, ie. the matching.
541 541

	
542 542
    /// @{
543 543

	
544 544
    ///\brief Returns the size of the current matching.
545 545
    ///
546 546
    ///Returns the size of the current matching. After \ref
547 547
    ///run() it returns the size of the maximum matching in the graph.
548 548
    int matchingSize() const {
549 549
      int size = 0;
550 550
      for (NodeIt n(_graph); n != INVALID; ++n) {
551 551
        if ((*_matching)[n] != INVALID) {
552 552
          ++size;
553 553
        }
554 554
      }
555 555
      return size / 2;
556 556
    }
557 557

	
558 558
    /// \brief Returns true when the edge is in the matching.
559 559
    ///
560 560
    /// Returns true when the edge is in the matching.
561 561
    bool matching(const Edge& edge) const {
562 562
      return edge == (*_matching)[_graph.u(edge)];
563 563
    }
... ...
@@ -1503,215 +1503,215 @@
1503 1503
        (*_blossom_data)[subblossoms[id]].next = next;
1504 1504
        (*_blossom_data)[subblossoms[id]].pred = pred;
1505 1505

	
1506 1506
      } else {
1507 1507

	
1508 1508
        for (int i = (ib + 1) % subblossoms.size();
1509 1509
             i != id; i = (i + 2) % subblossoms.size()) {
1510 1510
          int sb = subblossoms[i];
1511 1511
          int tb = subblossoms[(i + 1) % subblossoms.size()];
1512 1512
          (*_blossom_data)[sb].next =
1513 1513
            _graph.oppositeArc((*_blossom_data)[tb].next);
1514 1514
        }
1515 1515

	
1516 1516
        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1517 1517
          int sb = subblossoms[i];
1518 1518
          int tb = subblossoms[(i + 1) % subblossoms.size()];
1519 1519
          int ub = subblossoms[(i + 2) % subblossoms.size()];
1520 1520

	
1521 1521
          (*_blossom_data)[sb].status = ODD;
1522 1522
          matchedToOdd(sb);
1523 1523
          _tree_set->insert(sb, tree);
1524 1524
          (*_blossom_data)[sb].next = next;
1525 1525
          (*_blossom_data)[sb].pred =
1526 1526
            _graph.oppositeArc((*_blossom_data)[tb].next);
1527 1527

	
1528 1528
          (*_blossom_data)[tb].status = EVEN;
1529 1529
          matchedToEven(tb, tree);
1530 1530
          _tree_set->insert(tb, tree);
1531 1531
          (*_blossom_data)[tb].pred =
1532 1532
            (*_blossom_data)[tb].next =
1533 1533
            _graph.oppositeArc((*_blossom_data)[ub].next);
1534 1534
          next = (*_blossom_data)[ub].next;
1535 1535
        }
1536 1536

	
1537 1537
        (*_blossom_data)[subblossoms[ib]].status = ODD;
1538 1538
        matchedToOdd(subblossoms[ib]);
1539 1539
        _tree_set->insert(subblossoms[ib], tree);
1540 1540
        (*_blossom_data)[subblossoms[ib]].next = next;
1541 1541
        (*_blossom_data)[subblossoms[ib]].pred = pred;
1542 1542
      }
1543 1543
      _tree_set->erase(blossom);
1544 1544
    }
1545 1545

	
1546 1546
    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1547 1547
      if (_blossom_set->trivial(blossom)) {
1548 1548
        int bi = (*_node_index)[base];
1549 1549
        Value pot = (*_node_data)[bi].pot;
1550 1550

	
1551
        _matching->set(base, matching);
1551
        (*_matching)[base] = matching;
1552 1552
        _blossom_node_list.push_back(base);
1553
        _node_potential->set(base, pot);
1553
        (*_node_potential)[base] = pot;
1554 1554
      } else {
1555 1555

	
1556 1556
        Value pot = (*_blossom_data)[blossom].pot;
1557 1557
        int bn = _blossom_node_list.size();
1558 1558

	
1559 1559
        std::vector<int> subblossoms;
1560 1560
        _blossom_set->split(blossom, std::back_inserter(subblossoms));
1561 1561
        int b = _blossom_set->find(base);
1562 1562
        int ib = -1;
1563 1563
        for (int i = 0; i < int(subblossoms.size()); ++i) {
1564 1564
          if (subblossoms[i] == b) { ib = i; break; }
1565 1565
        }
1566 1566

	
1567 1567
        for (int i = 1; i < int(subblossoms.size()); i += 2) {
1568 1568
          int sb = subblossoms[(ib + i) % subblossoms.size()];
1569 1569
          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1570 1570

	
1571 1571
          Arc m = (*_blossom_data)[tb].next;
1572 1572
          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1573 1573
          extractBlossom(tb, _graph.source(m), m);
1574 1574
        }
1575 1575
        extractBlossom(subblossoms[ib], base, matching);
1576 1576

	
1577 1577
        int en = _blossom_node_list.size();
1578 1578

	
1579 1579
        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1580 1580
      }
1581 1581
    }
1582 1582

	
1583 1583
    void extractMatching() {
1584 1584
      std::vector<int> blossoms;
1585 1585
      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1586 1586
        blossoms.push_back(c);
1587 1587
      }
1588 1588

	
1589 1589
      for (int i = 0; i < int(blossoms.size()); ++i) {
1590 1590
        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
1591 1591

	
1592 1592
          Value offset = (*_blossom_data)[blossoms[i]].offset;
1593 1593
          (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1594 1594
          for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1595 1595
               n != INVALID; ++n) {
1596 1596
            (*_node_data)[(*_node_index)[n]].pot -= offset;
1597 1597
          }
1598 1598

	
1599 1599
          Arc matching = (*_blossom_data)[blossoms[i]].next;
1600 1600
          Node base = _graph.source(matching);
1601 1601
          extractBlossom(blossoms[i], base, matching);
1602 1602
        } else {
1603 1603
          Node base = (*_blossom_data)[blossoms[i]].base;
1604 1604
          extractBlossom(blossoms[i], base, INVALID);
1605 1605
        }
1606 1606
      }
1607 1607
    }
1608 1608

	
1609 1609
  public:
1610 1610

	
1611 1611
    /// \brief Constructor
1612 1612
    ///
1613 1613
    /// Constructor.
1614 1614
    MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1615 1615
      : _graph(graph), _weight(weight), _matching(0),
1616 1616
        _node_potential(0), _blossom_potential(), _blossom_node_list(),
1617 1617
        _node_num(0), _blossom_num(0),
1618 1618

	
1619 1619
        _blossom_index(0), _blossom_set(0), _blossom_data(0),
1620 1620
        _node_index(0), _node_heap_index(0), _node_data(0),
1621 1621
        _tree_set_index(0), _tree_set(0),
1622 1622

	
1623 1623
        _delta1_index(0), _delta1(0),
1624 1624
        _delta2_index(0), _delta2(0),
1625 1625
        _delta3_index(0), _delta3(0),
1626 1626
        _delta4_index(0), _delta4(0),
1627 1627

	
1628 1628
        _delta_sum() {}
1629 1629

	
1630 1630
    ~MaxWeightedMatching() {
1631 1631
      destroyStructures();
1632 1632
    }
1633 1633

	
1634 1634
    /// \name Execution control
1635 1635
    /// The simplest way to execute the algorithm is to use the
1636 1636
    /// \c run() member function.
1637 1637

	
1638 1638
    ///@{
1639 1639

	
1640 1640
    /// \brief Initialize the algorithm
1641 1641
    ///
1642 1642
    /// Initialize the algorithm
1643 1643
    void init() {
1644 1644
      createStructures();
1645 1645

	
1646 1646
      for (ArcIt e(_graph); e != INVALID; ++e) {
1647
        _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
1647
        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1648 1648
      }
1649 1649
      for (NodeIt n(_graph); n != INVALID; ++n) {
1650
        _delta1_index->set(n, _delta1->PRE_HEAP);
1650
        (*_delta1_index)[n] = _delta1->PRE_HEAP;
1651 1651
      }
1652 1652
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1653
        _delta3_index->set(e, _delta3->PRE_HEAP);
1653
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1654 1654
      }
1655 1655
      for (int i = 0; i < _blossom_num; ++i) {
1656
        _delta2_index->set(i, _delta2->PRE_HEAP);
1657
        _delta4_index->set(i, _delta4->PRE_HEAP);
1656
        (*_delta2_index)[i] = _delta2->PRE_HEAP;
1657
        (*_delta4_index)[i] = _delta4->PRE_HEAP;
1658 1658
      }
1659 1659

	
1660 1660
      int index = 0;
1661 1661
      for (NodeIt n(_graph); n != INVALID; ++n) {
1662 1662
        Value max = 0;
1663 1663
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1664 1664
          if (_graph.target(e) == n) continue;
1665 1665
          if ((dualScale * _weight[e]) / 2 > max) {
1666 1666
            max = (dualScale * _weight[e]) / 2;
1667 1667
          }
1668 1668
        }
1669
        _node_index->set(n, index);
1669
        (*_node_index)[n] = index;
1670 1670
        (*_node_data)[index].pot = max;
1671 1671
        _delta1->push(n, max);
1672 1672
        int blossom =
1673 1673
          _blossom_set->insert(n, std::numeric_limits<Value>::max());
1674 1674

	
1675 1675
        _tree_set->insert(blossom);
1676 1676

	
1677 1677
        (*_blossom_data)[blossom].status = EVEN;
1678 1678
        (*_blossom_data)[blossom].pred = INVALID;
1679 1679
        (*_blossom_data)[blossom].next = INVALID;
1680 1680
        (*_blossom_data)[blossom].pot = 0;
1681 1681
        (*_blossom_data)[blossom].offset = 0;
1682 1682
        ++index;
1683 1683
      }
1684 1684
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1685 1685
        int si = (*_node_index)[_graph.u(e)];
1686 1686
        int ti = (*_node_index)[_graph.v(e)];
1687 1687
        if (_graph.u(e) != _graph.v(e)) {
1688 1688
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1689 1689
                            dualScale * _weight[e]) / 2);
1690 1690
        }
1691 1691
      }
1692 1692
    }
1693 1693

	
1694 1694
    /// \brief Starts the algorithm
1695 1695
    ///
1696 1696
    /// Starts the algorithm
1697 1697
    void start() {
1698 1698
      enum OpType {
1699 1699
        D1, D2, D3, D4
1700 1700
      };
1701 1701

	
1702 1702
      int unmatched = _node_num;
1703 1703
      while (unmatched > 0) {
1704 1704
        Value d1 = !_delta1->empty() ?
1705 1705
          _delta1->prio() : std::numeric_limits<Value>::max();
1706 1706

	
1707 1707
        Value d2 = !_delta2->empty() ?
1708 1708
          _delta2->prio() : std::numeric_limits<Value>::max();
1709 1709

	
1710 1710
        Value d3 = !_delta3->empty() ?
1711 1711
          _delta3->prio() : std::numeric_limits<Value>::max();
1712 1712

	
1713 1713
        Value d4 = !_delta4->empty() ?
1714 1714
          _delta4->prio() : std::numeric_limits<Value>::max();
1715 1715

	
1716 1716
        _delta_sum = d1; OpType ot = D1;
1717 1717
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
... ...
@@ -2696,206 +2696,206 @@
2696 2696
        (*_blossom_data)[subblossoms[id]].next = next;
2697 2697
        (*_blossom_data)[subblossoms[id]].pred = pred;
2698 2698

	
2699 2699
      } else {
2700 2700

	
2701 2701
        for (int i = (ib + 1) % subblossoms.size();
2702 2702
             i != id; i = (i + 2) % subblossoms.size()) {
2703 2703
          int sb = subblossoms[i];
2704 2704
          int tb = subblossoms[(i + 1) % subblossoms.size()];
2705 2705
          (*_blossom_data)[sb].next =
2706 2706
            _graph.oppositeArc((*_blossom_data)[tb].next);
2707 2707
        }
2708 2708

	
2709 2709
        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2710 2710
          int sb = subblossoms[i];
2711 2711
          int tb = subblossoms[(i + 1) % subblossoms.size()];
2712 2712
          int ub = subblossoms[(i + 2) % subblossoms.size()];
2713 2713

	
2714 2714
          (*_blossom_data)[sb].status = ODD;
2715 2715
          matchedToOdd(sb);
2716 2716
          _tree_set->insert(sb, tree);
2717 2717
          (*_blossom_data)[sb].next = next;
2718 2718
          (*_blossom_data)[sb].pred =
2719 2719
            _graph.oppositeArc((*_blossom_data)[tb].next);
2720 2720

	
2721 2721
          (*_blossom_data)[tb].status = EVEN;
2722 2722
          matchedToEven(tb, tree);
2723 2723
          _tree_set->insert(tb, tree);
2724 2724
          (*_blossom_data)[tb].pred =
2725 2725
            (*_blossom_data)[tb].next =
2726 2726
            _graph.oppositeArc((*_blossom_data)[ub].next);
2727 2727
          next = (*_blossom_data)[ub].next;
2728 2728
        }
2729 2729

	
2730 2730
        (*_blossom_data)[subblossoms[ib]].status = ODD;
2731 2731
        matchedToOdd(subblossoms[ib]);
2732 2732
        _tree_set->insert(subblossoms[ib], tree);
2733 2733
        (*_blossom_data)[subblossoms[ib]].next = next;
2734 2734
        (*_blossom_data)[subblossoms[ib]].pred = pred;
2735 2735
      }
2736 2736
      _tree_set->erase(blossom);
2737 2737
    }
2738 2738

	
2739 2739
    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2740 2740
      if (_blossom_set->trivial(blossom)) {
2741 2741
        int bi = (*_node_index)[base];
2742 2742
        Value pot = (*_node_data)[bi].pot;
2743 2743

	
2744
        _matching->set(base, matching);
2744
        (*_matching)[base] = matching;
2745 2745
        _blossom_node_list.push_back(base);
2746
        _node_potential->set(base, pot);
2746
        (*_node_potential)[base] = pot;
2747 2747
      } else {
2748 2748

	
2749 2749
        Value pot = (*_blossom_data)[blossom].pot;
2750 2750
        int bn = _blossom_node_list.size();
2751 2751

	
2752 2752
        std::vector<int> subblossoms;
2753 2753
        _blossom_set->split(blossom, std::back_inserter(subblossoms));
2754 2754
        int b = _blossom_set->find(base);
2755 2755
        int ib = -1;
2756 2756
        for (int i = 0; i < int(subblossoms.size()); ++i) {
2757 2757
          if (subblossoms[i] == b) { ib = i; break; }
2758 2758
        }
2759 2759

	
2760 2760
        for (int i = 1; i < int(subblossoms.size()); i += 2) {
2761 2761
          int sb = subblossoms[(ib + i) % subblossoms.size()];
2762 2762
          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2763 2763

	
2764 2764
          Arc m = (*_blossom_data)[tb].next;
2765 2765
          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2766 2766
          extractBlossom(tb, _graph.source(m), m);
2767 2767
        }
2768 2768
        extractBlossom(subblossoms[ib], base, matching);
2769 2769

	
2770 2770
        int en = _blossom_node_list.size();
2771 2771

	
2772 2772
        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2773 2773
      }
2774 2774
    }
2775 2775

	
2776 2776
    void extractMatching() {
2777 2777
      std::vector<int> blossoms;
2778 2778
      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2779 2779
        blossoms.push_back(c);
2780 2780
      }
2781 2781

	
2782 2782
      for (int i = 0; i < int(blossoms.size()); ++i) {
2783 2783

	
2784 2784
        Value offset = (*_blossom_data)[blossoms[i]].offset;
2785 2785
        (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2786 2786
        for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2787 2787
             n != INVALID; ++n) {
2788 2788
          (*_node_data)[(*_node_index)[n]].pot -= offset;
2789 2789
        }
2790 2790

	
2791 2791
        Arc matching = (*_blossom_data)[blossoms[i]].next;
2792 2792
        Node base = _graph.source(matching);
2793 2793
        extractBlossom(blossoms[i], base, matching);
2794 2794
      }
2795 2795
    }
2796 2796

	
2797 2797
  public:
2798 2798

	
2799 2799
    /// \brief Constructor
2800 2800
    ///
2801 2801
    /// Constructor.
2802 2802
    MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2803 2803
      : _graph(graph), _weight(weight), _matching(0),
2804 2804
        _node_potential(0), _blossom_potential(), _blossom_node_list(),
2805 2805
        _node_num(0), _blossom_num(0),
2806 2806

	
2807 2807
        _blossom_index(0), _blossom_set(0), _blossom_data(0),
2808 2808
        _node_index(0), _node_heap_index(0), _node_data(0),
2809 2809
        _tree_set_index(0), _tree_set(0),
2810 2810

	
2811 2811
        _delta2_index(0), _delta2(0),
2812 2812
        _delta3_index(0), _delta3(0),
2813 2813
        _delta4_index(0), _delta4(0),
2814 2814

	
2815 2815
        _delta_sum() {}
2816 2816

	
2817 2817
    ~MaxWeightedPerfectMatching() {
2818 2818
      destroyStructures();
2819 2819
    }
2820 2820

	
2821 2821
    /// \name Execution control
2822 2822
    /// The simplest way to execute the algorithm is to use the
2823 2823
    /// \c run() member function.
2824 2824

	
2825 2825
    ///@{
2826 2826

	
2827 2827
    /// \brief Initialize the algorithm
2828 2828
    ///
2829 2829
    /// Initialize the algorithm
2830 2830
    void init() {
2831 2831
      createStructures();
2832 2832

	
2833 2833
      for (ArcIt e(_graph); e != INVALID; ++e) {
2834
        _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
2834
        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
2835 2835
      }
2836 2836
      for (EdgeIt e(_graph); e != INVALID; ++e) {
2837
        _delta3_index->set(e, _delta3->PRE_HEAP);
2837
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
2838 2838
      }
2839 2839
      for (int i = 0; i < _blossom_num; ++i) {
2840
        _delta2_index->set(i, _delta2->PRE_HEAP);
2841
        _delta4_index->set(i, _delta4->PRE_HEAP);
2840
        (*_delta2_index)[i] = _delta2->PRE_HEAP;
2841
        (*_delta4_index)[i] = _delta4->PRE_HEAP;
2842 2842
      }
2843 2843

	
2844 2844
      int index = 0;
2845 2845
      for (NodeIt n(_graph); n != INVALID; ++n) {
2846 2846
        Value max = - std::numeric_limits<Value>::max();
2847 2847
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2848 2848
          if (_graph.target(e) == n) continue;
2849 2849
          if ((dualScale * _weight[e]) / 2 > max) {
2850 2850
            max = (dualScale * _weight[e]) / 2;
2851 2851
          }
2852 2852
        }
2853
        _node_index->set(n, index);
2853
        (*_node_index)[n] = index;
2854 2854
        (*_node_data)[index].pot = max;
2855 2855
        int blossom =
2856 2856
          _blossom_set->insert(n, std::numeric_limits<Value>::max());
2857 2857

	
2858 2858
        _tree_set->insert(blossom);
2859 2859

	
2860 2860
        (*_blossom_data)[blossom].status = EVEN;
2861 2861
        (*_blossom_data)[blossom].pred = INVALID;
2862 2862
        (*_blossom_data)[blossom].next = INVALID;
2863 2863
        (*_blossom_data)[blossom].pot = 0;
2864 2864
        (*_blossom_data)[blossom].offset = 0;
2865 2865
        ++index;
2866 2866
      }
2867 2867
      for (EdgeIt e(_graph); e != INVALID; ++e) {
2868 2868
        int si = (*_node_index)[_graph.u(e)];
2869 2869
        int ti = (*_node_index)[_graph.v(e)];
2870 2870
        if (_graph.u(e) != _graph.v(e)) {
2871 2871
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
2872 2872
                            dualScale * _weight[e]) / 2);
2873 2873
        }
2874 2874
      }
2875 2875
    }
2876 2876

	
2877 2877
    /// \brief Starts the algorithm
2878 2878
    ///
2879 2879
    /// Starts the algorithm
2880 2880
    bool start() {
2881 2881
      enum OpType {
2882 2882
        D2, D3, D4
2883 2883
      };
2884 2884

	
2885 2885
      int unmatched = _node_num;
2886 2886
      while (unmatched > 0) {
2887 2887
        Value d2 = !_delta2->empty() ?
2888 2888
          _delta2->prio() : std::numeric_limits<Value>::max();
2889 2889

	
2890 2890
        Value d3 = !_delta3->empty() ?
2891 2891
          _delta3->prio() : std::numeric_limits<Value>::max();
2892 2892

	
2893 2893
        Value d4 = !_delta4->empty() ?
2894 2894
          _delta4->prio() : std::numeric_limits<Value>::max();
2895 2895

	
2896 2896
        _delta_sum = d2; OpType ot = D2;
2897 2897
        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
2898 2898
        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
2899 2899

	
2900 2900
        if (_delta_sum == std::numeric_limits<Value>::max()) {
2901 2901
          return false;
Show white space 96 line context
... ...
@@ -248,168 +248,168 @@
248 248
      if (local_pred) {
249 249
        delete _pred;
250 250
      }
251 251
      if (_arc_order) {
252 252
        delete _arc_order;
253 253
      }
254 254
      if (_node_order) {
255 255
        delete _node_order;
256 256
      }
257 257
      if (_cost_arcs) {
258 258
        delete _cost_arcs;
259 259
      }
260 260
      if (_heap) {
261 261
        delete _heap;
262 262
      }
263 263
      if (_heap_cross_ref) {
264 264
        delete _heap_cross_ref;
265 265
      }
266 266
    }
267 267

	
268 268
    Arc prepare(Node node) {
269 269
      std::vector<Node> nodes;
270 270
      (*_node_order)[node] = _dual_node_list.size();
271 271
      StackLevel level;
272 272
      level.node_level = _dual_node_list.size();
273 273
      _dual_node_list.push_back(node);
274 274
      for (InArcIt it(*_digraph, node); it != INVALID; ++it) {
275 275
        Arc arc = it;
276 276
        Node source = _digraph->source(arc);
277 277
        Value value = (*_cost)[it];
278 278
        if (source == node || (*_node_order)[source] == -3) continue;
279 279
        if ((*_cost_arcs)[source].arc == INVALID) {
280 280
          (*_cost_arcs)[source].arc = arc;
281 281
          (*_cost_arcs)[source].value = value;
282 282
          nodes.push_back(source);
283 283
        } else {
284 284
          if ((*_cost_arcs)[source].value > value) {
285 285
            (*_cost_arcs)[source].arc = arc;
286 286
            (*_cost_arcs)[source].value = value;
287 287
          }
288 288
        }
289 289
      }
290 290
      CostArc minimum = (*_cost_arcs)[nodes[0]];
291 291
      for (int i = 1; i < int(nodes.size()); ++i) {
292 292
        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
293 293
          minimum = (*_cost_arcs)[nodes[i]];
294 294
        }
295 295
      }
296
      _arc_order->set(minimum.arc, _dual_variables.size());
296
      (*_arc_order)[minimum.arc] = _dual_variables.size();
297 297
      DualVariable var(_dual_node_list.size() - 1,
298 298
                       _dual_node_list.size(), minimum.value);
299 299
      _dual_variables.push_back(var);
300 300
      for (int i = 0; i < int(nodes.size()); ++i) {
301 301
        (*_cost_arcs)[nodes[i]].value -= minimum.value;
302 302
        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
303 303
        (*_cost_arcs)[nodes[i]].arc = INVALID;
304 304
      }
305 305
      level_stack.push_back(level);
306 306
      return minimum.arc;
307 307
    }
308 308

	
309 309
    Arc contract(Node node) {
310 310
      int node_bottom = bottom(node);
311 311
      std::vector<Node> nodes;
312 312
      while (!level_stack.empty() &&
313 313
             level_stack.back().node_level >= node_bottom) {
314 314
        for (int i = 0; i < int(level_stack.back().arcs.size()); ++i) {
315 315
          Arc arc = level_stack.back().arcs[i].arc;
316 316
          Node source = _digraph->source(arc);
317 317
          Value value = level_stack.back().arcs[i].value;
318 318
          if ((*_node_order)[source] >= node_bottom) continue;
319 319
          if ((*_cost_arcs)[source].arc == INVALID) {
320 320
            (*_cost_arcs)[source].arc = arc;
321 321
            (*_cost_arcs)[source].value = value;
322 322
            nodes.push_back(source);
323 323
          } else {
324 324
            if ((*_cost_arcs)[source].value > value) {
325 325
              (*_cost_arcs)[source].arc = arc;
326 326
              (*_cost_arcs)[source].value = value;
327 327
            }
328 328
          }
329 329
        }
330 330
        level_stack.pop_back();
331 331
      }
332 332
      CostArc minimum = (*_cost_arcs)[nodes[0]];
333 333
      for (int i = 1; i < int(nodes.size()); ++i) {
334 334
        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
335 335
          minimum = (*_cost_arcs)[nodes[i]];
336 336
        }
337 337
      }
338
      _arc_order->set(minimum.arc, _dual_variables.size());
338
      (*_arc_order)[minimum.arc] = _dual_variables.size();
339 339
      DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
340 340
      _dual_variables.push_back(var);
341 341
      StackLevel level;
342 342
      level.node_level = node_bottom;
343 343
      for (int i = 0; i < int(nodes.size()); ++i) {
344 344
        (*_cost_arcs)[nodes[i]].value -= minimum.value;
345 345
        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
346 346
        (*_cost_arcs)[nodes[i]].arc = INVALID;
347 347
      }
348 348
      level_stack.push_back(level);
349 349
      return minimum.arc;
350 350
    }
351 351

	
352 352
    int bottom(Node node) {
353 353
      int k = level_stack.size() - 1;
354 354
      while (level_stack[k].node_level > (*_node_order)[node]) {
355 355
        --k;
356 356
      }
357 357
      return level_stack[k].node_level;
358 358
    }
359 359

	
360 360
    void finalize(Arc arc) {
361 361
      Node node = _digraph->target(arc);
362 362
      _heap->push(node, (*_arc_order)[arc]);
363 363
      _pred->set(node, arc);
364 364
      while (!_heap->empty()) {
365 365
        Node source = _heap->top();
366 366
        _heap->pop();
367
        _node_order->set(source, -1);
367
        (*_node_order)[source] = -1;
368 368
        for (OutArcIt it(*_digraph, source); it != INVALID; ++it) {
369 369
          if ((*_arc_order)[it] < 0) continue;
370 370
          Node target = _digraph->target(it);
371 371
          switch(_heap->state(target)) {
372 372
          case Heap::PRE_HEAP:
373 373
            _heap->push(target, (*_arc_order)[it]);
374 374
            _pred->set(target, it);
375 375
            break;
376 376
          case Heap::IN_HEAP:
377 377
            if ((*_arc_order)[it] < (*_heap)[target]) {
378 378
              _heap->decrease(target, (*_arc_order)[it]);
379 379
              _pred->set(target, it);
380 380
            }
381 381
            break;
382 382
          case Heap::POST_HEAP:
383 383
            break;
384 384
          }
385 385
        }
386 386
        _arborescence->set((*_pred)[source], true);
387 387
      }
388 388
    }
389 389

	
390 390

	
391 391
  public:
392 392

	
393 393
    /// \name Named template parameters
394 394

	
395 395
    /// @{
396 396

	
397 397
    template <class T>
398 398
    struct DefArborescenceMapTraits : public Traits {
399 399
      typedef T ArborescenceMap;
400 400
      static ArborescenceMap *createArborescenceMap(const Digraph &)
401 401
      {
402 402
        LEMON_ASSERT(false, "ArborescenceMap is not initialized");
403 403
        return 0; // ignore warnings
404 404
      }
405 405
    };
406 406

	
407 407
    /// \brief \ref named-templ-param "Named parameter" for
408 408
    /// setting ArborescenceMap type
409 409
    ///
410 410
    /// \ref named-templ-param "Named parameter" for setting
411 411
    /// ArborescenceMap type
412 412
    template <class T>
413 413
    struct DefArborescenceMap
414 414
      : public MinCostArborescence<Digraph, CostMap,
415 415
                                   DefArborescenceMapTraits<T> > {
... ...
@@ -605,103 +605,103 @@
605 605
      ///
606 606
      /// Increment operator.
607 607
      DualIt& operator++() {
608 608
        ++_index;
609 609
        return *this;
610 610
      }
611 611

	
612 612
      /// \brief Validity checking
613 613
      ///
614 614
      /// Checks whether the iterator is invalid.
615 615
      bool operator==(Invalid) const {
616 616
        return _index == _last;
617 617
      }
618 618

	
619 619
      /// \brief Validity checking
620 620
      ///
621 621
      /// Checks whether the iterator is valid.
622 622
      bool operator!=(Invalid) const {
623 623
        return _index != _last;
624 624
      }
625 625

	
626 626
    private:
627 627
      const MinCostArborescence* _algorithm;
628 628
      int _index, _last;
629 629
    };
630 630

	
631 631
    /// @}
632 632

	
633 633
    /// \name Execution control
634 634
    /// The simplest way to execute the algorithm is to use
635 635
    /// one of the member functions called \c run(...). \n
636 636
    /// If you need more control on the execution,
637 637
    /// first you must call \ref init(), then you can add several
638 638
    /// source nodes with \ref addSource().
639 639
    /// Finally \ref start() will perform the arborescence
640 640
    /// computation.
641 641

	
642 642
    ///@{
643 643

	
644 644
    /// \brief Initializes the internal data structures.
645 645
    ///
646 646
    /// Initializes the internal data structures.
647 647
    ///
648 648
    void init() {
649 649
      createStructures();
650 650
      _heap->clear();
651 651
      for (NodeIt it(*_digraph); it != INVALID; ++it) {
652 652
        (*_cost_arcs)[it].arc = INVALID;
653
        _node_order->set(it, -3);
654
        _heap_cross_ref->set(it, Heap::PRE_HEAP);
653
        (*_node_order)[it] = -3;
654
        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
655 655
        _pred->set(it, INVALID);
656 656
      }
657 657
      for (ArcIt it(*_digraph); it != INVALID; ++it) {
658 658
        _arborescence->set(it, false);
659
        _arc_order->set(it, -1);
659
        (*_arc_order)[it] = -1;
660 660
      }
661 661
      _dual_node_list.clear();
662 662
      _dual_variables.clear();
663 663
    }
664 664

	
665 665
    /// \brief Adds a new source node.
666 666
    ///
667 667
    /// Adds a new source node to the algorithm.
668 668
    void addSource(Node source) {
669 669
      std::vector<Node> nodes;
670 670
      nodes.push_back(source);
671 671
      while (!nodes.empty()) {
672 672
        Node node = nodes.back();
673 673
        nodes.pop_back();
674 674
        for (OutArcIt it(*_digraph, node); it != INVALID; ++it) {
675 675
          Node target = _digraph->target(it);
676 676
          if ((*_node_order)[target] == -3) {
677 677
            (*_node_order)[target] = -2;
678 678
            nodes.push_back(target);
679 679
            queue.push_back(target);
680 680
          }
681 681
        }
682 682
      }
683 683
      (*_node_order)[source] = -1;
684 684
    }
685 685

	
686 686
    /// \brief Processes the next node in the priority queue.
687 687
    ///
688 688
    /// Processes the next node in the priority queue.
689 689
    ///
690 690
    /// \return The processed node.
691 691
    ///
692 692
    /// \warning The queue must not be empty!
693 693
    Node processNextNode() {
694 694
      Node node = queue.back();
695 695
      queue.pop_back();
696 696
      if ((*_node_order)[node] == -2) {
697 697
        Arc arc = prepare(node);
698 698
        Node source = _digraph->source(arc);
699 699
        while ((*_node_order)[source] != -1) {
700 700
          if ((*_node_order)[source] >= 0) {
701 701
            arc = contract(source);
702 702
          } else {
703 703
            arc = prepare(source);
704 704
          }
705 705
          source = _digraph->source(arc);
706 706
        }
707 707
        finalize(arc);
Show white space 96 line context
... ...
@@ -359,523 +359,523 @@
359 359
      }
360 360
      _level = &elevator;
361 361
      return *this;
362 362
    }
363 363

	
364 364
    /// \brief Returns a const reference to the elevator.
365 365
    ///
366 366
    /// Returns a const reference to the elevator.
367 367
    ///
368 368
    /// \pre Either \ref run() or \ref init() must be called before
369 369
    /// using this function.
370 370
    const Elevator& elevator() const {
371 371
      return *_level;
372 372
    }
373 373

	
374 374
    /// \brief Sets the tolerance used by algorithm.
375 375
    ///
376 376
    /// Sets the tolerance used by algorithm.
377 377
    Preflow& tolerance(const Tolerance& tolerance) const {
378 378
      _tolerance = tolerance;
379 379
      return *this;
380 380
    }
381 381

	
382 382
    /// \brief Returns a const reference to the tolerance.
383 383
    ///
384 384
    /// Returns a const reference to the tolerance.
385 385
    const Tolerance& tolerance() const {
386 386
      return tolerance;
387 387
    }
388 388

	
389 389
    /// \name Execution Control
390 390
    /// The simplest way to execute the preflow algorithm is to use
391 391
    /// \ref run() or \ref runMinCut().\n
392 392
    /// If you need more control on the initial solution or the execution,
393 393
    /// first you have to call one of the \ref init() functions, then
394 394
    /// \ref startFirstPhase() and if you need it \ref startSecondPhase().
395 395

	
396 396
    ///@{
397 397

	
398 398
    /// \brief Initializes the internal data structures.
399 399
    ///
400 400
    /// Initializes the internal data structures and sets the initial
401 401
    /// flow to zero on each arc.
402 402
    void init() {
403 403
      createStructures();
404 404

	
405 405
      _phase = true;
406 406
      for (NodeIt n(_graph); n != INVALID; ++n) {
407
        _excess->set(n, 0);
407
        (*_excess)[n] = 0;
408 408
      }
409 409

	
410 410
      for (ArcIt e(_graph); e != INVALID; ++e) {
411 411
        _flow->set(e, 0);
412 412
      }
413 413

	
414 414
      typename Digraph::template NodeMap<bool> reached(_graph, false);
415 415

	
416 416
      _level->initStart();
417 417
      _level->initAddItem(_target);
418 418

	
419 419
      std::vector<Node> queue;
420
      reached.set(_source, true);
420
      reached[_source] = true;
421 421

	
422 422
      queue.push_back(_target);
423
      reached.set(_target, true);
423
      reached[_target] = true;
424 424
      while (!queue.empty()) {
425 425
        _level->initNewLevel();
426 426
        std::vector<Node> nqueue;
427 427
        for (int i = 0; i < int(queue.size()); ++i) {
428 428
          Node n = queue[i];
429 429
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
430 430
            Node u = _graph.source(e);
431 431
            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
432
              reached.set(u, true);
432
              reached[u] = true;
433 433
              _level->initAddItem(u);
434 434
              nqueue.push_back(u);
435 435
            }
436 436
          }
437 437
        }
438 438
        queue.swap(nqueue);
439 439
      }
440 440
      _level->initFinish();
441 441

	
442 442
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
443 443
        if (_tolerance.positive((*_capacity)[e])) {
444 444
          Node u = _graph.target(e);
445 445
          if ((*_level)[u] == _level->maxLevel()) continue;
446 446
          _flow->set(e, (*_capacity)[e]);
447
          _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
447
          (*_excess)[u] += (*_capacity)[e];
448 448
          if (u != _target && !_level->active(u)) {
449 449
            _level->activate(u);
450 450
          }
451 451
        }
452 452
      }
453 453
    }
454 454

	
455 455
    /// \brief Initializes the internal data structures using the
456 456
    /// given flow map.
457 457
    ///
458 458
    /// Initializes the internal data structures and sets the initial
459 459
    /// flow to the given \c flowMap. The \c flowMap should contain a
460 460
    /// flow or at least a preflow, i.e. at each node excluding the
461 461
    /// source node the incoming flow should greater or equal to the
462 462
    /// outgoing flow.
463 463
    /// \return \c false if the given \c flowMap is not a preflow.
464 464
    template <typename FlowMap>
465 465
    bool init(const FlowMap& flowMap) {
466 466
      createStructures();
467 467

	
468 468
      for (ArcIt e(_graph); e != INVALID; ++e) {
469 469
        _flow->set(e, flowMap[e]);
470 470
      }
471 471

	
472 472
      for (NodeIt n(_graph); n != INVALID; ++n) {
473 473
        Value excess = 0;
474 474
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
475 475
          excess += (*_flow)[e];
476 476
        }
477 477
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
478 478
          excess -= (*_flow)[e];
479 479
        }
480 480
        if (excess < 0 && n != _source) return false;
481
        _excess->set(n, excess);
481
        (*_excess)[n] = excess;
482 482
      }
483 483

	
484 484
      typename Digraph::template NodeMap<bool> reached(_graph, false);
485 485

	
486 486
      _level->initStart();
487 487
      _level->initAddItem(_target);
488 488

	
489 489
      std::vector<Node> queue;
490
      reached.set(_source, true);
490
      reached[_source] = true;
491 491

	
492 492
      queue.push_back(_target);
493
      reached.set(_target, true);
493
      reached[_target] = true;
494 494
      while (!queue.empty()) {
495 495
        _level->initNewLevel();
496 496
        std::vector<Node> nqueue;
497 497
        for (int i = 0; i < int(queue.size()); ++i) {
498 498
          Node n = queue[i];
499 499
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
500 500
            Node u = _graph.source(e);
501 501
            if (!reached[u] &&
502 502
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
503
              reached.set(u, true);
503
              reached[u] = true;
504 504
              _level->initAddItem(u);
505 505
              nqueue.push_back(u);
506 506
            }
507 507
          }
508 508
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
509 509
            Node v = _graph.target(e);
510 510
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
511
              reached.set(v, true);
511
              reached[v] = true;
512 512
              _level->initAddItem(v);
513 513
              nqueue.push_back(v);
514 514
            }
515 515
          }
516 516
        }
517 517
        queue.swap(nqueue);
518 518
      }
519 519
      _level->initFinish();
520 520

	
521 521
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
522 522
        Value rem = (*_capacity)[e] - (*_flow)[e];
523 523
        if (_tolerance.positive(rem)) {
524 524
          Node u = _graph.target(e);
525 525
          if ((*_level)[u] == _level->maxLevel()) continue;
526 526
          _flow->set(e, (*_capacity)[e]);
527
          _excess->set(u, (*_excess)[u] + rem);
527
          (*_excess)[u] += rem;
528 528
          if (u != _target && !_level->active(u)) {
529 529
            _level->activate(u);
530 530
          }
531 531
        }
532 532
      }
533 533
      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
534 534
        Value rem = (*_flow)[e];
535 535
        if (_tolerance.positive(rem)) {
536 536
          Node v = _graph.source(e);
537 537
          if ((*_level)[v] == _level->maxLevel()) continue;
538 538
          _flow->set(e, 0);
539
          _excess->set(v, (*_excess)[v] + rem);
539
          (*_excess)[v] += rem;
540 540
          if (v != _target && !_level->active(v)) {
541 541
            _level->activate(v);
542 542
          }
543 543
        }
544 544
      }
545 545
      return true;
546 546
    }
547 547

	
548 548
    /// \brief Starts the first phase of the preflow algorithm.
549 549
    ///
550 550
    /// The preflow algorithm consists of two phases, this method runs
551 551
    /// the first phase. After the first phase the maximum flow value
552 552
    /// and a minimum value cut can already be computed, although a
553 553
    /// maximum flow is not yet obtained. So after calling this method
554 554
    /// \ref flowValue() returns the value of a maximum flow and \ref
555 555
    /// minCut() returns a minimum cut.
556 556
    /// \pre One of the \ref init() functions must be called before
557 557
    /// using this function.
558 558
    void startFirstPhase() {
559 559
      _phase = true;
560 560

	
561 561
      Node n = _level->highestActive();
562 562
      int level = _level->highestActiveLevel();
563 563
      while (n != INVALID) {
564 564
        int num = _node_num;
565 565

	
566 566
        while (num > 0 && n != INVALID) {
567 567
          Value excess = (*_excess)[n];
568 568
          int new_level = _level->maxLevel();
569 569

	
570 570
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
571 571
            Value rem = (*_capacity)[e] - (*_flow)[e];
572 572
            if (!_tolerance.positive(rem)) continue;
573 573
            Node v = _graph.target(e);
574 574
            if ((*_level)[v] < level) {
575 575
              if (!_level->active(v) && v != _target) {
576 576
                _level->activate(v);
577 577
              }
578 578
              if (!_tolerance.less(rem, excess)) {
579 579
                _flow->set(e, (*_flow)[e] + excess);
580
                _excess->set(v, (*_excess)[v] + excess);
580
                (*_excess)[v] += excess;
581 581
                excess = 0;
582 582
                goto no_more_push_1;
583 583
              } else {
584 584
                excess -= rem;
585
                _excess->set(v, (*_excess)[v] + rem);
585
                (*_excess)[v] += rem;
586 586
                _flow->set(e, (*_capacity)[e]);
587 587
              }
588 588
            } else if (new_level > (*_level)[v]) {
589 589
              new_level = (*_level)[v];
590 590
            }
591 591
          }
592 592

	
593 593
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
594 594
            Value rem = (*_flow)[e];
595 595
            if (!_tolerance.positive(rem)) continue;
596 596
            Node v = _graph.source(e);
597 597
            if ((*_level)[v] < level) {
598 598
              if (!_level->active(v) && v != _target) {
599 599
                _level->activate(v);
600 600
              }
601 601
              if (!_tolerance.less(rem, excess)) {
602 602
                _flow->set(e, (*_flow)[e] - excess);
603
                _excess->set(v, (*_excess)[v] + excess);
603
                (*_excess)[v] += excess;
604 604
                excess = 0;
605 605
                goto no_more_push_1;
606 606
              } else {
607 607
                excess -= rem;
608
                _excess->set(v, (*_excess)[v] + rem);
608
                (*_excess)[v] += rem;
609 609
                _flow->set(e, 0);
610 610
              }
611 611
            } else if (new_level > (*_level)[v]) {
612 612
              new_level = (*_level)[v];
613 613
            }
614 614
          }
615 615

	
616 616
        no_more_push_1:
617 617

	
618
          _excess->set(n, excess);
618
          (*_excess)[n] = excess;
619 619

	
620 620
          if (excess != 0) {
621 621
            if (new_level + 1 < _level->maxLevel()) {
622 622
              _level->liftHighestActive(new_level + 1);
623 623
            } else {
624 624
              _level->liftHighestActiveToTop();
625 625
            }
626 626
            if (_level->emptyLevel(level)) {
627 627
              _level->liftToTop(level);
628 628
            }
629 629
          } else {
630 630
            _level->deactivate(n);
631 631
          }
632 632

	
633 633
          n = _level->highestActive();
634 634
          level = _level->highestActiveLevel();
635 635
          --num;
636 636
        }
637 637

	
638 638
        num = _node_num * 20;
639 639
        while (num > 0 && n != INVALID) {
640 640
          Value excess = (*_excess)[n];
641 641
          int new_level = _level->maxLevel();
642 642

	
643 643
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
644 644
            Value rem = (*_capacity)[e] - (*_flow)[e];
645 645
            if (!_tolerance.positive(rem)) continue;
646 646
            Node v = _graph.target(e);
647 647
            if ((*_level)[v] < level) {
648 648
              if (!_level->active(v) && v != _target) {
649 649
                _level->activate(v);
650 650
              }
651 651
              if (!_tolerance.less(rem, excess)) {
652 652
                _flow->set(e, (*_flow)[e] + excess);
653
                _excess->set(v, (*_excess)[v] + excess);
653
                (*_excess)[v] += excess;
654 654
                excess = 0;
655 655
                goto no_more_push_2;
656 656
              } else {
657 657
                excess -= rem;
658
                _excess->set(v, (*_excess)[v] + rem);
658
                (*_excess)[v] += rem;
659 659
                _flow->set(e, (*_capacity)[e]);
660 660
              }
661 661
            } else if (new_level > (*_level)[v]) {
662 662
              new_level = (*_level)[v];
663 663
            }
664 664
          }
665 665

	
666 666
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
667 667
            Value rem = (*_flow)[e];
668 668
            if (!_tolerance.positive(rem)) continue;
669 669
            Node v = _graph.source(e);
670 670
            if ((*_level)[v] < level) {
671 671
              if (!_level->active(v) && v != _target) {
672 672
                _level->activate(v);
673 673
              }
674 674
              if (!_tolerance.less(rem, excess)) {
675 675
                _flow->set(e, (*_flow)[e] - excess);
676
                _excess->set(v, (*_excess)[v] + excess);
676
                (*_excess)[v] += excess;
677 677
                excess = 0;
678 678
                goto no_more_push_2;
679 679
              } else {
680 680
                excess -= rem;
681
                _excess->set(v, (*_excess)[v] + rem);
681
                (*_excess)[v] += rem;
682 682
                _flow->set(e, 0);
683 683
              }
684 684
            } else if (new_level > (*_level)[v]) {
685 685
              new_level = (*_level)[v];
686 686
            }
687 687
          }
688 688

	
689 689
        no_more_push_2:
690 690

	
691
          _excess->set(n, excess);
691
          (*_excess)[n] = excess;
692 692

	
693 693
          if (excess != 0) {
694 694
            if (new_level + 1 < _level->maxLevel()) {
695 695
              _level->liftActiveOn(level, new_level + 1);
696 696
            } else {
697 697
              _level->liftActiveToTop(level);
698 698
            }
699 699
            if (_level->emptyLevel(level)) {
700 700
              _level->liftToTop(level);
701 701
            }
702 702
          } else {
703 703
            _level->deactivate(n);
704 704
          }
705 705

	
706 706
          while (level >= 0 && _level->activeFree(level)) {
707 707
            --level;
708 708
          }
709 709
          if (level == -1) {
710 710
            n = _level->highestActive();
711 711
            level = _level->highestActiveLevel();
712 712
          } else {
713 713
            n = _level->activeOn(level);
714 714
          }
715 715
          --num;
716 716
        }
717 717
      }
718 718
    }
719 719

	
720 720
    /// \brief Starts the second phase of the preflow algorithm.
721 721
    ///
722 722
    /// The preflow algorithm consists of two phases, this method runs
723 723
    /// the second phase. After calling one of the \ref init() functions
724 724
    /// and \ref startFirstPhase() and then \ref startSecondPhase(),
725 725
    /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
726 726
    /// value of a maximum flow, \ref minCut() returns a minimum cut
727 727
    /// \pre One of the \ref init() functions and \ref startFirstPhase()
728 728
    /// must be called before using this function.
729 729
    void startSecondPhase() {
730 730
      _phase = false;
731 731

	
732 732
      typename Digraph::template NodeMap<bool> reached(_graph);
733 733
      for (NodeIt n(_graph); n != INVALID; ++n) {
734
        reached.set(n, (*_level)[n] < _level->maxLevel());
734
        reached[n] = (*_level)[n] < _level->maxLevel();
735 735
      }
736 736

	
737 737
      _level->initStart();
738 738
      _level->initAddItem(_source);
739 739

	
740 740
      std::vector<Node> queue;
741 741
      queue.push_back(_source);
742
      reached.set(_source, true);
742
      reached[_source] = true;
743 743

	
744 744
      while (!queue.empty()) {
745 745
        _level->initNewLevel();
746 746
        std::vector<Node> nqueue;
747 747
        for (int i = 0; i < int(queue.size()); ++i) {
748 748
          Node n = queue[i];
749 749
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
750 750
            Node v = _graph.target(e);
751 751
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
752
              reached.set(v, true);
752
              reached[v] = true;
753 753
              _level->initAddItem(v);
754 754
              nqueue.push_back(v);
755 755
            }
756 756
          }
757 757
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
758 758
            Node u = _graph.source(e);
759 759
            if (!reached[u] &&
760 760
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
761
              reached.set(u, true);
761
              reached[u] = true;
762 762
              _level->initAddItem(u);
763 763
              nqueue.push_back(u);
764 764
            }
765 765
          }
766 766
        }
767 767
        queue.swap(nqueue);
768 768
      }
769 769
      _level->initFinish();
770 770

	
771 771
      for (NodeIt n(_graph); n != INVALID; ++n) {
772 772
        if (!reached[n]) {
773 773
          _level->dirtyTopButOne(n);
774 774
        } else if ((*_excess)[n] > 0 && _target != n) {
775 775
          _level->activate(n);
776 776
        }
777 777
      }
778 778

	
779 779
      Node n;
780 780
      while ((n = _level->highestActive()) != INVALID) {
781 781
        Value excess = (*_excess)[n];
782 782
        int level = _level->highestActiveLevel();
783 783
        int new_level = _level->maxLevel();
784 784

	
785 785
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
786 786
          Value rem = (*_capacity)[e] - (*_flow)[e];
787 787
          if (!_tolerance.positive(rem)) continue;
788 788
          Node v = _graph.target(e);
789 789
          if ((*_level)[v] < level) {
790 790
            if (!_level->active(v) && v != _source) {
791 791
              _level->activate(v);
792 792
            }
793 793
            if (!_tolerance.less(rem, excess)) {
794 794
              _flow->set(e, (*_flow)[e] + excess);
795
              _excess->set(v, (*_excess)[v] + excess);
795
              (*_excess)[v] += excess;
796 796
              excess = 0;
797 797
              goto no_more_push;
798 798
            } else {
799 799
              excess -= rem;
800
              _excess->set(v, (*_excess)[v] + rem);
800
              (*_excess)[v] += rem;
801 801
              _flow->set(e, (*_capacity)[e]);
802 802
            }
803 803
          } else if (new_level > (*_level)[v]) {
804 804
            new_level = (*_level)[v];
805 805
          }
806 806
        }
807 807

	
808 808
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
809 809
          Value rem = (*_flow)[e];
810 810
          if (!_tolerance.positive(rem)) continue;
811 811
          Node v = _graph.source(e);
812 812
          if ((*_level)[v] < level) {
813 813
            if (!_level->active(v) && v != _source) {
814 814
              _level->activate(v);
815 815
            }
816 816
            if (!_tolerance.less(rem, excess)) {
817 817
              _flow->set(e, (*_flow)[e] - excess);
818
              _excess->set(v, (*_excess)[v] + excess);
818
              (*_excess)[v] += excess;
819 819
              excess = 0;
820 820
              goto no_more_push;
821 821
            } else {
822 822
              excess -= rem;
823
              _excess->set(v, (*_excess)[v] + rem);
823
              (*_excess)[v] += rem;
824 824
              _flow->set(e, 0);
825 825
            }
826 826
          } else if (new_level > (*_level)[v]) {
827 827
            new_level = (*_level)[v];
828 828
          }
829 829
        }
830 830

	
831 831
      no_more_push:
832 832

	
833
        _excess->set(n, excess);
833
        (*_excess)[n] = excess;
834 834

	
835 835
        if (excess != 0) {
836 836
          if (new_level + 1 < _level->maxLevel()) {
837 837
            _level->liftHighestActive(new_level + 1);
838 838
          } else {
839 839
            // Calculation error
840 840
            _level->liftHighestActiveToTop();
841 841
          }
842 842
          if (_level->emptyLevel(level)) {
843 843
            // Calculation error
844 844
            _level->liftToTop(level);
845 845
          }
846 846
        } else {
847 847
          _level->deactivate(n);
848 848
        }
849 849

	
850 850
      }
851 851
    }
852 852

	
853 853
    /// \brief Runs the preflow algorithm.
854 854
    ///
855 855
    /// Runs the preflow algorithm.
856 856
    /// \note pf.run() is just a shortcut of the following code.
857 857
    /// \code
858 858
    ///   pf.init();
859 859
    ///   pf.startFirstPhase();
860 860
    ///   pf.startSecondPhase();
861 861
    /// \endcode
862 862
    void run() {
863 863
      init();
864 864
      startFirstPhase();
865 865
      startSecondPhase();
866 866
    }
867 867

	
868 868
    /// \brief Runs the preflow algorithm to compute the minimum cut.
869 869
    ///
870 870
    /// Runs the preflow algorithm to compute the minimum cut.
871 871
    /// \note pf.runMinCut() is just a shortcut of the following code.
872 872
    /// \code
873 873
    ///   pf.init();
874 874
    ///   pf.startFirstPhase();
875 875
    /// \endcode
876 876
    void runMinCut() {
877 877
      init();
878 878
      startFirstPhase();
879 879
    }
880 880

	
881 881
    /// @}
Show white space 96 line context
... ...
@@ -54,94 +54,94 @@
54 54
  std::vector<concepts::Digraph::Arc> ws;
55 55
  std::vector<concepts::Graph::Edge> uws;
56 56

	
57 57
  kruskal(g, r, ws.begin());
58 58
  kruskal(ug, ur, uws.begin());
59 59
}
60 60

	
61 61
int main() {
62 62

	
63 63
  typedef ListGraph::Node Node;
64 64
  typedef ListGraph::Edge Edge;
65 65
  typedef ListGraph::NodeIt NodeIt;
66 66
  typedef ListGraph::ArcIt ArcIt;
67 67

	
68 68
  ListGraph G;
69 69

	
70 70
  Node s=G.addNode();
71 71
  Node v1=G.addNode();
72 72
  Node v2=G.addNode();
73 73
  Node v3=G.addNode();
74 74
  Node v4=G.addNode();
75 75
  Node t=G.addNode();
76 76

	
77 77
  Edge e1 = G.addEdge(s, v1);
78 78
  Edge e2 = G.addEdge(s, v2);
79 79
  Edge e3 = G.addEdge(v1, v2);
80 80
  Edge e4 = G.addEdge(v2, v1);
81 81
  Edge e5 = G.addEdge(v1, v3);
82 82
  Edge e6 = G.addEdge(v3, v2);
83 83
  Edge e7 = G.addEdge(v2, v4);
84 84
  Edge e8 = G.addEdge(v4, v3);
85 85
  Edge e9 = G.addEdge(v3, t);
86 86
  Edge e10 = G.addEdge(v4, t);
87 87

	
88 88
  typedef ListGraph::EdgeMap<int> ECostMap;
89 89
  typedef ListGraph::EdgeMap<bool> EBoolMap;
90 90

	
91 91
  ECostMap edge_cost_map(G, 2);
92 92
  EBoolMap tree_map(G);
93 93

	
94 94

	
95 95
  //Test with const map.
96 96
  check(kruskal(G, ConstMap<ListGraph::Edge,int>(2), tree_map)==10,
97 97
        "Total cost should be 10");
98 98
  //Test with an edge map (filled with uniform costs).
99 99
  check(kruskal(G, edge_cost_map, tree_map)==10,
100 100
        "Total cost should be 10");
101 101

	
102
  edge_cost_map.set(e1, -10);
103
  edge_cost_map.set(e2, -9);
104
  edge_cost_map.set(e3, -8);
105
  edge_cost_map.set(e4, -7);
106
  edge_cost_map.set(e5, -6);
107
  edge_cost_map.set(e6, -5);
108
  edge_cost_map.set(e7, -4);
109
  edge_cost_map.set(e8, -3);
110
  edge_cost_map.set(e9, -2);
111
  edge_cost_map.set(e10, -1);
102
  edge_cost_map[e1] = -10;
103
  edge_cost_map[e2] = -9;
104
  edge_cost_map[e3] = -8;
105
  edge_cost_map[e4] = -7;
106
  edge_cost_map[e5] = -6;
107
  edge_cost_map[e6] = -5;
108
  edge_cost_map[e7] = -4;
109
  edge_cost_map[e8] = -3;
110
  edge_cost_map[e9] = -2;
111
  edge_cost_map[e10] = -1;
112 112

	
113 113
  vector<Edge> tree_edge_vec(5);
114 114

	
115 115
  //Test with a edge map and inserter.
116 116
  check(kruskal(G, edge_cost_map,
117 117
                 tree_edge_vec.begin())
118 118
        ==-31,
119 119
        "Total cost should be -31.");
120 120

	
121 121
  tree_edge_vec.clear();
122 122

	
123 123
  check(kruskal(G, edge_cost_map,
124 124
                back_inserter(tree_edge_vec))
125 125
        ==-31,
126 126
        "Total cost should be -31.");
127 127

	
128 128
//   tree_edge_vec.clear();
129 129

	
130 130
//   //The above test could also be coded like this:
131 131
//   check(kruskal(G,
132 132
//                 makeKruskalMapInput(G, edge_cost_map),
133 133
//                 makeKruskalSequenceOutput(back_inserter(tree_edge_vec)))
134 134
//         ==-31,
135 135
//         "Total cost should be -31.");
136 136

	
137 137
  check(tree_edge_vec.size()==5,"The tree should have 5 edges.");
138 138

	
139 139
  check(tree_edge_vec[0]==e1 &&
140 140
        tree_edge_vec[1]==e2 &&
141 141
        tree_edge_vec[2]==e5 &&
142 142
        tree_edge_vec[3]==e7 &&
143 143
        tree_edge_vec[4]==e9,
144 144
        "Wrong tree.");
145 145

	
146 146
  return 0;
147 147
}
0 comments (0 inline)