gravatar
deba@inf.elte.hu
deba@inf.elte.hu
Fix critical bug in preflow (#372) The wrong transition between the bound decrease and highest active heuristics caused the bug. The last node chosen in bound decrease mode is used in the first iteration in highest active mode.
0 1 0
default
1 file changed with 24 insertions and 20 deletions:
↑ Collapse diff ↑
Ignore white space 192 line context
... ...
@@ -465,349 +465,353 @@
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 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 490
      reached[_source] = true;
491 491

	
492 492
      queue.push_back(_target);
493 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 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 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 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 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
      Node n = _level->highestActive();
562
      int level = _level->highestActiveLevel();
563
      while (n != INVALID) {
561
      while (true) {
564 562
        int num = _node_num;
565 563

	
566
        while (num > 0 && n != INVALID) {
564
        Node n = INVALID;
565
        int level = -1;
566

	
567
        while (num > 0) {
568
          n = _level->highestActive();
569
          if (n == INVALID) goto first_phase_done;
570
          level = _level->highestActiveLevel();
571
          --num;
572
          
567 573
          Value excess = (*_excess)[n];
568 574
          int new_level = _level->maxLevel();
569 575

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

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

	
616 622
        no_more_push_1:
617 623

	
618 624
          (*_excess)[n] = excess;
619 625

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

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

	
638 640
        num = _node_num * 20;
639
        while (num > 0 && n != INVALID) {
641
        while (num > 0) {
642
          while (level >= 0 && _level->activeFree(level)) {
643
            --level;
644
          }
645
          if (level == -1) {
646
            n = _level->highestActive();
647
            level = _level->highestActiveLevel();
648
            if (n == INVALID) goto first_phase_done;
649
          } else {
650
            n = _level->activeOn(level);
651
          }
652
          --num;
653

	
640 654
          Value excess = (*_excess)[n];
641 655
          int new_level = _level->maxLevel();
642 656

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

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

	
689 703
        no_more_push_2:
690 704

	
691 705
          (*_excess)[n] = excess;
692 706

	
693 707
          if (excess != 0) {
694 708
            if (new_level + 1 < _level->maxLevel()) {
695 709
              _level->liftActiveOn(level, new_level + 1);
696 710
            } else {
697 711
              _level->liftActiveToTop(level);
698 712
            }
699 713
            if (_level->emptyLevel(level)) {
700 714
              _level->liftToTop(level);
701 715
            }
702 716
          } else {
703 717
            _level->deactivate(n);
704 718
          }
705

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

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

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

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

	
740 744
      std::vector<Node> queue;
741 745
      queue.push_back(_source);
742 746
      reached[_source] = true;
743 747

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

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

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

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

	
808 812
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
809 813
          Value rem = (*_flow)[e];
810 814
          if (!_tolerance.positive(rem)) continue;
811 815
          Node v = _graph.source(e);
812 816
          if ((*_level)[v] < level) {
813 817
            if (!_level->active(v) && v != _source) {
0 comments (0 inline)