gravatar
alpar (Alpar Juttner)
alpar@cs.elte.hu
Merge #356
0 2 0
merge default
2 files changed with 71 insertions and 2 deletions:
↑ Collapse diff ↑
Ignore white space 6 line context
... ...
@@ -429,2949 +429,3009 @@
429 429
    ///@{
430 430

	
431 431
    /// \brief Set the initial matching to the empty matching.
432 432
    ///
433 433
    /// This function sets the initial matching to the empty matching.
434 434
    void init() {
435 435
      createStructures();
436 436
      for(NodeIt n(_graph); n != INVALID; ++n) {
437 437
        (*_matching)[n] = INVALID;
438 438
        (*_status)[n] = UNMATCHED;
439 439
      }
440 440
    }
441 441

	
442 442
    /// \brief Find an initial matching in a greedy way.
443 443
    ///
444 444
    /// This function finds an initial matching in a greedy way.
445 445
    void greedyInit() {
446 446
      createStructures();
447 447
      for (NodeIt n(_graph); n != INVALID; ++n) {
448 448
        (*_matching)[n] = INVALID;
449 449
        (*_status)[n] = UNMATCHED;
450 450
      }
451 451
      for (NodeIt n(_graph); n != INVALID; ++n) {
452 452
        if ((*_matching)[n] == INVALID) {
453 453
          for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
454 454
            Node v = _graph.target(a);
455 455
            if ((*_matching)[v] == INVALID && v != n) {
456 456
              (*_matching)[n] = a;
457 457
              (*_status)[n] = MATCHED;
458 458
              (*_matching)[v] = _graph.oppositeArc(a);
459 459
              (*_status)[v] = MATCHED;
460 460
              break;
461 461
            }
462 462
          }
463 463
        }
464 464
      }
465 465
    }
466 466

	
467 467

	
468 468
    /// \brief Initialize the matching from a map.
469 469
    ///
470 470
    /// This function initializes the matching from a \c bool valued edge
471 471
    /// map. This map should have the property that there are no two incident
472 472
    /// edges with \c true value, i.e. it really contains a matching.
473 473
    /// \return \c true if the map contains a matching.
474 474
    template <typename MatchingMap>
475 475
    bool matchingInit(const MatchingMap& matching) {
476 476
      createStructures();
477 477

	
478 478
      for (NodeIt n(_graph); n != INVALID; ++n) {
479 479
        (*_matching)[n] = INVALID;
480 480
        (*_status)[n] = UNMATCHED;
481 481
      }
482 482
      for(EdgeIt e(_graph); e!=INVALID; ++e) {
483 483
        if (matching[e]) {
484 484

	
485 485
          Node u = _graph.u(e);
486 486
          if ((*_matching)[u] != INVALID) return false;
487 487
          (*_matching)[u] = _graph.direct(e, true);
488 488
          (*_status)[u] = MATCHED;
489 489

	
490 490
          Node v = _graph.v(e);
491 491
          if ((*_matching)[v] != INVALID) return false;
492 492
          (*_matching)[v] = _graph.direct(e, false);
493 493
          (*_status)[v] = MATCHED;
494 494
        }
495 495
      }
496 496
      return true;
497 497
    }
498 498

	
499 499
    /// \brief Start Edmonds' algorithm
500 500
    ///
501 501
    /// This function runs the original Edmonds' algorithm.
502 502
    ///
503 503
    /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
504 504
    /// called before using this function.
505 505
    void startSparse() {
506 506
      for(NodeIt n(_graph); n != INVALID; ++n) {
507 507
        if ((*_status)[n] == UNMATCHED) {
508 508
          (*_blossom_rep)[_blossom_set->insert(n)] = n;
509 509
          _tree_set->insert(n);
510 510
          (*_status)[n] = EVEN;
511 511
          processSparse(n);
512 512
        }
513 513
      }
514 514
    }
515 515

	
516 516
    /// \brief Start Edmonds' algorithm with a heuristic improvement
517 517
    /// for dense graphs
518 518
    ///
519 519
    /// This function runs Edmonds' algorithm with a heuristic of postponing
520 520
    /// shrinks, therefore resulting in a faster algorithm for dense graphs.
521 521
    ///
522 522
    /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
523 523
    /// called before using this function.
524 524
    void startDense() {
525 525
      for(NodeIt n(_graph); n != INVALID; ++n) {
526 526
        if ((*_status)[n] == UNMATCHED) {
527 527
          (*_blossom_rep)[_blossom_set->insert(n)] = n;
528 528
          _tree_set->insert(n);
529 529
          (*_status)[n] = EVEN;
530 530
          processDense(n);
531 531
        }
532 532
      }
533 533
    }
534 534

	
535 535

	
536 536
    /// \brief Run Edmonds' algorithm
537 537
    ///
538 538
    /// This function runs Edmonds' algorithm. An additional heuristic of
539 539
    /// postponing shrinks is used for relatively dense graphs
540 540
    /// (for which <tt>m>=2*n</tt> holds).
541 541
    void run() {
542 542
      if (countEdges(_graph) < 2 * countNodes(_graph)) {
543 543
        greedyInit();
544 544
        startSparse();
545 545
      } else {
546 546
        init();
547 547
        startDense();
548 548
      }
549 549
    }
550 550

	
551 551
    /// @}
552 552

	
553 553
    /// \name Primal Solution
554 554
    /// Functions to get the primal solution, i.e. the maximum matching.
555 555

	
556 556
    /// @{
557 557

	
558 558
    /// \brief Return the size (cardinality) of the matching.
559 559
    ///
560 560
    /// This function returns the size (cardinality) of the current matching.
561 561
    /// After run() it returns the size of the maximum matching in the graph.
562 562
    int matchingSize() const {
563 563
      int size = 0;
564 564
      for (NodeIt n(_graph); n != INVALID; ++n) {
565 565
        if ((*_matching)[n] != INVALID) {
566 566
          ++size;
567 567
        }
568 568
      }
569 569
      return size / 2;
570 570
    }
571 571

	
572 572
    /// \brief Return \c true if the given edge is in the matching.
573 573
    ///
574 574
    /// This function returns \c true if the given edge is in the current
575 575
    /// matching.
576 576
    bool matching(const Edge& edge) const {
577 577
      return edge == (*_matching)[_graph.u(edge)];
578 578
    }
579 579

	
580 580
    /// \brief Return the matching arc (or edge) incident to the given node.
581 581
    ///
582 582
    /// This function returns the matching arc (or edge) incident to the
583 583
    /// given node in the current matching or \c INVALID if the node is
584 584
    /// not covered by the matching.
585 585
    Arc matching(const Node& n) const {
586 586
      return (*_matching)[n];
587 587
    }
588 588

	
589 589
    /// \brief Return a const reference to the matching map.
590 590
    ///
591 591
    /// This function returns a const reference to a node map that stores
592 592
    /// the matching arc (or edge) incident to each node.
593 593
    const MatchingMap& matchingMap() const {
594 594
      return *_matching;
595 595
    }
596 596

	
597 597
    /// \brief Return the mate of the given node.
598 598
    ///
599 599
    /// This function returns the mate of the given node in the current
600 600
    /// matching or \c INVALID if the node is not covered by the matching.
601 601
    Node mate(const Node& n) const {
602 602
      return (*_matching)[n] != INVALID ?
603 603
        _graph.target((*_matching)[n]) : INVALID;
604 604
    }
605 605

	
606 606
    /// @}
607 607

	
608 608
    /// \name Dual Solution
609 609
    /// Functions to get the dual solution, i.e. the Gallai-Edmonds
610 610
    /// decomposition.
611 611

	
612 612
    /// @{
613 613

	
614 614
    /// \brief Return the status of the given node in the Edmonds-Gallai
615 615
    /// decomposition.
616 616
    ///
617 617
    /// This function returns the \ref Status "status" of the given node
618 618
    /// in the Edmonds-Gallai decomposition.
619 619
    Status status(const Node& n) const {
620 620
      return (*_status)[n];
621 621
    }
622 622

	
623 623
    /// \brief Return a const reference to the status map, which stores
624 624
    /// the Edmonds-Gallai decomposition.
625 625
    ///
626 626
    /// This function returns a const reference to a node map that stores the
627 627
    /// \ref Status "status" of each node in the Edmonds-Gallai decomposition.
628 628
    const StatusMap& statusMap() const {
629 629
      return *_status;
630 630
    }
631 631

	
632 632
    /// \brief Return \c true if the given node is in the barrier.
633 633
    ///
634 634
    /// This function returns \c true if the given node is in the barrier.
635 635
    bool barrier(const Node& n) const {
636 636
      return (*_status)[n] == ODD;
637 637
    }
638 638

	
639 639
    /// @}
640 640

	
641 641
  };
642 642

	
643 643
  /// \ingroup matching
644 644
  ///
645 645
  /// \brief Weighted matching in general graphs
646 646
  ///
647 647
  /// This class provides an efficient implementation of Edmond's
648 648
  /// maximum weighted matching algorithm. The implementation is based
649 649
  /// on extensive use of priority queues and provides
650 650
  /// \f$O(nm\log n)\f$ time complexity.
651 651
  ///
652 652
  /// The maximum weighted matching problem is to find a subset of the
653 653
  /// edges in an undirected graph with maximum overall weight for which
654 654
  /// each node has at most one incident edge.
655 655
  /// It can be formulated with the following linear program.
656 656
  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
657 657
  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
658 658
      \quad \forall B\in\mathcal{O}\f] */
659 659
  /// \f[x_e \ge 0\quad \forall e\in E\f]
660 660
  /// \f[\max \sum_{e\in E}x_ew_e\f]
661 661
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
662 662
  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
663 663
  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
664 664
  /// subsets of the nodes.
665 665
  ///
666 666
  /// The algorithm calculates an optimal matching and a proof of the
667 667
  /// optimality. The solution of the dual problem can be used to check
668 668
  /// the result of the algorithm. The dual linear problem is the
669 669
  /// following.
670 670
  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
671 671
      z_B \ge w_{uv} \quad \forall uv\in E\f] */
672 672
  /// \f[y_u \ge 0 \quad \forall u \in V\f]
673 673
  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
674 674
  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
675 675
      \frac{\vert B \vert - 1}{2}z_B\f] */
676 676
  ///
677 677
  /// The algorithm can be executed with the run() function.
678 678
  /// After it the matching (the primal solution) and the dual solution
679 679
  /// can be obtained using the query functions and the
680 680
  /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
681 681
  /// which is able to iterate on the nodes of a blossom.
682 682
  /// If the value type is integer, then the dual solution is multiplied
683 683
  /// by \ref MaxWeightedMatching::dualScale "4".
684 684
  ///
685 685
  /// \tparam GR The undirected graph type the algorithm runs on.
686 686
  /// \tparam WM The type edge weight map. The default type is
687 687
  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
688 688
#ifdef DOXYGEN
689 689
  template <typename GR, typename WM>
690 690
#else
691 691
  template <typename GR,
692 692
            typename WM = typename GR::template EdgeMap<int> >
693 693
#endif
694 694
  class MaxWeightedMatching {
695 695
  public:
696 696

	
697 697
    /// The graph type of the algorithm
698 698
    typedef GR Graph;
699 699
    /// The type of the edge weight map
700 700
    typedef WM WeightMap;
701 701
    /// The value type of the edge weights
702 702
    typedef typename WeightMap::Value Value;
703 703

	
704 704
    /// The type of the matching map
705 705
    typedef typename Graph::template NodeMap<typename Graph::Arc>
706 706
    MatchingMap;
707 707

	
708 708
    /// \brief Scaling factor for dual solution
709 709
    ///
710 710
    /// Scaling factor for dual solution. It is equal to 4 or 1
711 711
    /// according to the value type.
712 712
    static const int dualScale =
713 713
      std::numeric_limits<Value>::is_integer ? 4 : 1;
714 714

	
715 715
  private:
716 716

	
717 717
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
718 718

	
719 719
    typedef typename Graph::template NodeMap<Value> NodePotential;
720 720
    typedef std::vector<Node> BlossomNodeList;
721 721

	
722 722
    struct BlossomVariable {
723 723
      int begin, end;
724 724
      Value value;
725 725

	
726 726
      BlossomVariable(int _begin, int _end, Value _value)
727 727
        : begin(_begin), end(_end), value(_value) {}
728 728

	
729 729
    };
730 730

	
731 731
    typedef std::vector<BlossomVariable> BlossomPotential;
732 732

	
733 733
    const Graph& _graph;
734 734
    const WeightMap& _weight;
735 735

	
736 736
    MatchingMap* _matching;
737 737

	
738 738
    NodePotential* _node_potential;
739 739

	
740 740
    BlossomPotential _blossom_potential;
741 741
    BlossomNodeList _blossom_node_list;
742 742

	
743 743
    int _node_num;
744 744
    int _blossom_num;
745 745

	
746 746
    typedef RangeMap<int> IntIntMap;
747 747

	
748 748
    enum Status {
749 749
      EVEN = -1, MATCHED = 0, ODD = 1
750 750
    };
751 751

	
752 752
    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
753 753
    struct BlossomData {
754 754
      int tree;
755 755
      Status status;
756 756
      Arc pred, next;
757 757
      Value pot, offset;
758 758
      Node base;
759 759
    };
760 760

	
761 761
    IntNodeMap *_blossom_index;
762 762
    BlossomSet *_blossom_set;
763 763
    RangeMap<BlossomData>* _blossom_data;
764 764

	
765 765
    IntNodeMap *_node_index;
766 766
    IntArcMap *_node_heap_index;
767 767

	
768 768
    struct NodeData {
769 769

	
770 770
      NodeData(IntArcMap& node_heap_index)
771 771
        : heap(node_heap_index) {}
772 772

	
773 773
      int blossom;
774 774
      Value pot;
775 775
      BinHeap<Value, IntArcMap> heap;
776 776
      std::map<int, Arc> heap_index;
777 777

	
778 778
      int tree;
779 779
    };
780 780

	
781 781
    RangeMap<NodeData>* _node_data;
782 782

	
783 783
    typedef ExtendFindEnum<IntIntMap> TreeSet;
784 784

	
785 785
    IntIntMap *_tree_set_index;
786 786
    TreeSet *_tree_set;
787 787

	
788 788
    IntNodeMap *_delta1_index;
789 789
    BinHeap<Value, IntNodeMap> *_delta1;
790 790

	
791 791
    IntIntMap *_delta2_index;
792 792
    BinHeap<Value, IntIntMap> *_delta2;
793 793

	
794 794
    IntEdgeMap *_delta3_index;
795 795
    BinHeap<Value, IntEdgeMap> *_delta3;
796 796

	
797 797
    IntIntMap *_delta4_index;
798 798
    BinHeap<Value, IntIntMap> *_delta4;
799 799

	
800 800
    Value _delta_sum;
801 801
    int _unmatched;
802 802

	
803 803
    typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
804 804
    FractionalMatching *_fractional;
805 805

	
806 806
    void createStructures() {
807 807
      _node_num = countNodes(_graph);
808 808
      _blossom_num = _node_num * 3 / 2;
809 809

	
810 810
      if (!_matching) {
811 811
        _matching = new MatchingMap(_graph);
812 812
      }
813

	
813 814
      if (!_node_potential) {
814 815
        _node_potential = new NodePotential(_graph);
815 816
      }
817

	
816 818
      if (!_blossom_set) {
817 819
        _blossom_index = new IntNodeMap(_graph);
818 820
        _blossom_set = new BlossomSet(*_blossom_index);
819 821
        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
822
      } else if (_blossom_data->size() != _blossom_num) {
823
        delete _blossom_data;
824
        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
820 825
      }
821 826

	
822 827
      if (!_node_index) {
823 828
        _node_index = new IntNodeMap(_graph);
824 829
        _node_heap_index = new IntArcMap(_graph);
825 830
        _node_data = new RangeMap<NodeData>(_node_num,
826
                                              NodeData(*_node_heap_index));
831
                                            NodeData(*_node_heap_index));
832
      } else {
833
        delete _node_data;
834
        _node_data = new RangeMap<NodeData>(_node_num,
835
                                            NodeData(*_node_heap_index));
827 836
      }
828 837

	
829 838
      if (!_tree_set) {
830 839
        _tree_set_index = new IntIntMap(_blossom_num);
831 840
        _tree_set = new TreeSet(*_tree_set_index);
841
      } else {
842
        _tree_set_index->resize(_blossom_num);
832 843
      }
844

	
833 845
      if (!_delta1) {
834 846
        _delta1_index = new IntNodeMap(_graph);
835 847
        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
836 848
      }
849

	
837 850
      if (!_delta2) {
838 851
        _delta2_index = new IntIntMap(_blossom_num);
839 852
        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
853
      } else {
854
        _delta2_index->resize(_blossom_num);
840 855
      }
856

	
841 857
      if (!_delta3) {
842 858
        _delta3_index = new IntEdgeMap(_graph);
843 859
        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
844 860
      }
861

	
845 862
      if (!_delta4) {
846 863
        _delta4_index = new IntIntMap(_blossom_num);
847 864
        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
865
      } else {
866
        _delta4_index->resize(_blossom_num);
848 867
      }
849 868
    }
850 869

	
851 870
    void destroyStructures() {
852 871
      if (_matching) {
853 872
        delete _matching;
854 873
      }
855 874
      if (_node_potential) {
856 875
        delete _node_potential;
857 876
      }
858 877
      if (_blossom_set) {
859 878
        delete _blossom_index;
860 879
        delete _blossom_set;
861 880
        delete _blossom_data;
862 881
      }
863 882

	
864 883
      if (_node_index) {
865 884
        delete _node_index;
866 885
        delete _node_heap_index;
867 886
        delete _node_data;
868 887
      }
869 888

	
870 889
      if (_tree_set) {
871 890
        delete _tree_set_index;
872 891
        delete _tree_set;
873 892
      }
874 893
      if (_delta1) {
875 894
        delete _delta1_index;
876 895
        delete _delta1;
877 896
      }
878 897
      if (_delta2) {
879 898
        delete _delta2_index;
880 899
        delete _delta2;
881 900
      }
882 901
      if (_delta3) {
883 902
        delete _delta3_index;
884 903
        delete _delta3;
885 904
      }
886 905
      if (_delta4) {
887 906
        delete _delta4_index;
888 907
        delete _delta4;
889 908
      }
890 909
    }
891 910

	
892 911
    void matchedToEven(int blossom, int tree) {
893 912
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
894 913
        _delta2->erase(blossom);
895 914
      }
896 915

	
897 916
      if (!_blossom_set->trivial(blossom)) {
898 917
        (*_blossom_data)[blossom].pot -=
899 918
          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
900 919
      }
901 920

	
902 921
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
903 922
           n != INVALID; ++n) {
904 923

	
905 924
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
906 925
        int ni = (*_node_index)[n];
907 926

	
908 927
        (*_node_data)[ni].heap.clear();
909 928
        (*_node_data)[ni].heap_index.clear();
910 929

	
911 930
        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
912 931

	
913 932
        _delta1->push(n, (*_node_data)[ni].pot);
914 933

	
915 934
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
916 935
          Node v = _graph.source(e);
917 936
          int vb = _blossom_set->find(v);
918 937
          int vi = (*_node_index)[v];
919 938

	
920 939
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
921 940
            dualScale * _weight[e];
922 941

	
923 942
          if ((*_blossom_data)[vb].status == EVEN) {
924 943
            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
925 944
              _delta3->push(e, rw / 2);
926 945
            }
927 946
          } else {
928 947
            typename std::map<int, Arc>::iterator it =
929 948
              (*_node_data)[vi].heap_index.find(tree);
930 949

	
931 950
            if (it != (*_node_data)[vi].heap_index.end()) {
932 951
              if ((*_node_data)[vi].heap[it->second] > rw) {
933 952
                (*_node_data)[vi].heap.replace(it->second, e);
934 953
                (*_node_data)[vi].heap.decrease(e, rw);
935 954
                it->second = e;
936 955
              }
937 956
            } else {
938 957
              (*_node_data)[vi].heap.push(e, rw);
939 958
              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
940 959
            }
941 960

	
942 961
            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
943 962
              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
944 963

	
945 964
              if ((*_blossom_data)[vb].status == MATCHED) {
946 965
                if (_delta2->state(vb) != _delta2->IN_HEAP) {
947 966
                  _delta2->push(vb, _blossom_set->classPrio(vb) -
948 967
                               (*_blossom_data)[vb].offset);
949 968
                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
950 969
                           (*_blossom_data)[vb].offset) {
951 970
                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
952 971
                                   (*_blossom_data)[vb].offset);
953 972
                }
954 973
              }
955 974
            }
956 975
          }
957 976
        }
958 977
      }
959 978
      (*_blossom_data)[blossom].offset = 0;
960 979
    }
961 980

	
962 981
    void matchedToOdd(int blossom) {
963 982
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
964 983
        _delta2->erase(blossom);
965 984
      }
966 985
      (*_blossom_data)[blossom].offset += _delta_sum;
967 986
      if (!_blossom_set->trivial(blossom)) {
968 987
        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
969 988
                      (*_blossom_data)[blossom].offset);
970 989
      }
971 990
    }
972 991

	
973 992
    void evenToMatched(int blossom, int tree) {
974 993
      if (!_blossom_set->trivial(blossom)) {
975 994
        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
976 995
      }
977 996

	
978 997
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
979 998
           n != INVALID; ++n) {
980 999
        int ni = (*_node_index)[n];
981 1000
        (*_node_data)[ni].pot -= _delta_sum;
982 1001

	
983 1002
        _delta1->erase(n);
984 1003

	
985 1004
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
986 1005
          Node v = _graph.source(e);
987 1006
          int vb = _blossom_set->find(v);
988 1007
          int vi = (*_node_index)[v];
989 1008

	
990 1009
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
991 1010
            dualScale * _weight[e];
992 1011

	
993 1012
          if (vb == blossom) {
994 1013
            if (_delta3->state(e) == _delta3->IN_HEAP) {
995 1014
              _delta3->erase(e);
996 1015
            }
997 1016
          } else if ((*_blossom_data)[vb].status == EVEN) {
998 1017

	
999 1018
            if (_delta3->state(e) == _delta3->IN_HEAP) {
1000 1019
              _delta3->erase(e);
1001 1020
            }
1002 1021

	
1003 1022
            int vt = _tree_set->find(vb);
1004 1023

	
1005 1024
            if (vt != tree) {
1006 1025

	
1007 1026
              Arc r = _graph.oppositeArc(e);
1008 1027

	
1009 1028
              typename std::map<int, Arc>::iterator it =
1010 1029
                (*_node_data)[ni].heap_index.find(vt);
1011 1030

	
1012 1031
              if (it != (*_node_data)[ni].heap_index.end()) {
1013 1032
                if ((*_node_data)[ni].heap[it->second] > rw) {
1014 1033
                  (*_node_data)[ni].heap.replace(it->second, r);
1015 1034
                  (*_node_data)[ni].heap.decrease(r, rw);
1016 1035
                  it->second = r;
1017 1036
                }
1018 1037
              } else {
1019 1038
                (*_node_data)[ni].heap.push(r, rw);
1020 1039
                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1021 1040
              }
1022 1041

	
1023 1042
              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1024 1043
                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1025 1044

	
1026 1045
                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1027 1046
                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1028 1047
                               (*_blossom_data)[blossom].offset);
1029 1048
                } else if ((*_delta2)[blossom] >
1030 1049
                           _blossom_set->classPrio(blossom) -
1031 1050
                           (*_blossom_data)[blossom].offset){
1032 1051
                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1033 1052
                                   (*_blossom_data)[blossom].offset);
1034 1053
                }
1035 1054
              }
1036 1055
            }
1037 1056
          } else {
1038 1057

	
1039 1058
            typename std::map<int, Arc>::iterator it =
1040 1059
              (*_node_data)[vi].heap_index.find(tree);
1041 1060

	
1042 1061
            if (it != (*_node_data)[vi].heap_index.end()) {
1043 1062
              (*_node_data)[vi].heap.erase(it->second);
1044 1063
              (*_node_data)[vi].heap_index.erase(it);
1045 1064
              if ((*_node_data)[vi].heap.empty()) {
1046 1065
                _blossom_set->increase(v, std::numeric_limits<Value>::max());
1047 1066
              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1048 1067
                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1049 1068
              }
1050 1069

	
1051 1070
              if ((*_blossom_data)[vb].status == MATCHED) {
1052 1071
                if (_blossom_set->classPrio(vb) ==
1053 1072
                    std::numeric_limits<Value>::max()) {
1054 1073
                  _delta2->erase(vb);
1055 1074
                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1056 1075
                           (*_blossom_data)[vb].offset) {
1057 1076
                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
1058 1077
                                   (*_blossom_data)[vb].offset);
1059 1078
                }
1060 1079
              }
1061 1080
            }
1062 1081
          }
1063 1082
        }
1064 1083
      }
1065 1084
    }
1066 1085

	
1067 1086
    void oddToMatched(int blossom) {
1068 1087
      (*_blossom_data)[blossom].offset -= _delta_sum;
1069 1088

	
1070 1089
      if (_blossom_set->classPrio(blossom) !=
1071 1090
          std::numeric_limits<Value>::max()) {
1072 1091
        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1073 1092
                      (*_blossom_data)[blossom].offset);
1074 1093
      }
1075 1094

	
1076 1095
      if (!_blossom_set->trivial(blossom)) {
1077 1096
        _delta4->erase(blossom);
1078 1097
      }
1079 1098
    }
1080 1099

	
1081 1100
    void oddToEven(int blossom, int tree) {
1082 1101
      if (!_blossom_set->trivial(blossom)) {
1083 1102
        _delta4->erase(blossom);
1084 1103
        (*_blossom_data)[blossom].pot -=
1085 1104
          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1086 1105
      }
1087 1106

	
1088 1107
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1089 1108
           n != INVALID; ++n) {
1090 1109
        int ni = (*_node_index)[n];
1091 1110

	
1092 1111
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
1093 1112

	
1094 1113
        (*_node_data)[ni].heap.clear();
1095 1114
        (*_node_data)[ni].heap_index.clear();
1096 1115
        (*_node_data)[ni].pot +=
1097 1116
          2 * _delta_sum - (*_blossom_data)[blossom].offset;
1098 1117

	
1099 1118
        _delta1->push(n, (*_node_data)[ni].pot);
1100 1119

	
1101 1120
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
1102 1121
          Node v = _graph.source(e);
1103 1122
          int vb = _blossom_set->find(v);
1104 1123
          int vi = (*_node_index)[v];
1105 1124

	
1106 1125
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1107 1126
            dualScale * _weight[e];
1108 1127

	
1109 1128
          if ((*_blossom_data)[vb].status == EVEN) {
1110 1129
            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1111 1130
              _delta3->push(e, rw / 2);
1112 1131
            }
1113 1132
          } else {
1114 1133

	
1115 1134
            typename std::map<int, Arc>::iterator it =
1116 1135
              (*_node_data)[vi].heap_index.find(tree);
1117 1136

	
1118 1137
            if (it != (*_node_data)[vi].heap_index.end()) {
1119 1138
              if ((*_node_data)[vi].heap[it->second] > rw) {
1120 1139
                (*_node_data)[vi].heap.replace(it->second, e);
1121 1140
                (*_node_data)[vi].heap.decrease(e, rw);
1122 1141
                it->second = e;
1123 1142
              }
1124 1143
            } else {
1125 1144
              (*_node_data)[vi].heap.push(e, rw);
1126 1145
              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1127 1146
            }
1128 1147

	
1129 1148
            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1130 1149
              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1131 1150

	
1132 1151
              if ((*_blossom_data)[vb].status == MATCHED) {
1133 1152
                if (_delta2->state(vb) != _delta2->IN_HEAP) {
1134 1153
                  _delta2->push(vb, _blossom_set->classPrio(vb) -
1135 1154
                               (*_blossom_data)[vb].offset);
1136 1155
                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1137 1156
                           (*_blossom_data)[vb].offset) {
1138 1157
                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1139 1158
                                   (*_blossom_data)[vb].offset);
1140 1159
                }
1141 1160
              }
1142 1161
            }
1143 1162
          }
1144 1163
        }
1145 1164
      }
1146 1165
      (*_blossom_data)[blossom].offset = 0;
1147 1166
    }
1148 1167

	
1149 1168
    void alternatePath(int even, int tree) {
1150 1169
      int odd;
1151 1170

	
1152 1171
      evenToMatched(even, tree);
1153 1172
      (*_blossom_data)[even].status = MATCHED;
1154 1173

	
1155 1174
      while ((*_blossom_data)[even].pred != INVALID) {
1156 1175
        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
1157 1176
        (*_blossom_data)[odd].status = MATCHED;
1158 1177
        oddToMatched(odd);
1159 1178
        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1160 1179

	
1161 1180
        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
1162 1181
        (*_blossom_data)[even].status = MATCHED;
1163 1182
        evenToMatched(even, tree);
1164 1183
        (*_blossom_data)[even].next =
1165 1184
          _graph.oppositeArc((*_blossom_data)[odd].pred);
1166 1185
      }
1167 1186

	
1168 1187
    }
1169 1188

	
1170 1189
    void destroyTree(int tree) {
1171 1190
      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1172 1191
        if ((*_blossom_data)[b].status == EVEN) {
1173 1192
          (*_blossom_data)[b].status = MATCHED;
1174 1193
          evenToMatched(b, tree);
1175 1194
        } else if ((*_blossom_data)[b].status == ODD) {
1176 1195
          (*_blossom_data)[b].status = MATCHED;
1177 1196
          oddToMatched(b);
1178 1197
        }
1179 1198
      }
1180 1199
      _tree_set->eraseClass(tree);
1181 1200
    }
1182 1201

	
1183 1202

	
1184 1203
    void unmatchNode(const Node& node) {
1185 1204
      int blossom = _blossom_set->find(node);
1186 1205
      int tree = _tree_set->find(blossom);
1187 1206

	
1188 1207
      alternatePath(blossom, tree);
1189 1208
      destroyTree(tree);
1190 1209

	
1191 1210
      (*_blossom_data)[blossom].base = node;
1192 1211
      (*_blossom_data)[blossom].next = INVALID;
1193 1212
    }
1194 1213

	
1195 1214
    void augmentOnEdge(const Edge& edge) {
1196 1215

	
1197 1216
      int left = _blossom_set->find(_graph.u(edge));
1198 1217
      int right = _blossom_set->find(_graph.v(edge));
1199 1218

	
1200 1219
      int left_tree = _tree_set->find(left);
1201 1220
      alternatePath(left, left_tree);
1202 1221
      destroyTree(left_tree);
1203 1222

	
1204 1223
      int right_tree = _tree_set->find(right);
1205 1224
      alternatePath(right, right_tree);
1206 1225
      destroyTree(right_tree);
1207 1226

	
1208 1227
      (*_blossom_data)[left].next = _graph.direct(edge, true);
1209 1228
      (*_blossom_data)[right].next = _graph.direct(edge, false);
1210 1229
    }
1211 1230

	
1212 1231
    void augmentOnArc(const Arc& arc) {
1213 1232

	
1214 1233
      int left = _blossom_set->find(_graph.source(arc));
1215 1234
      int right = _blossom_set->find(_graph.target(arc));
1216 1235

	
1217 1236
      (*_blossom_data)[left].status = MATCHED;
1218 1237

	
1219 1238
      int right_tree = _tree_set->find(right);
1220 1239
      alternatePath(right, right_tree);
1221 1240
      destroyTree(right_tree);
1222 1241

	
1223 1242
      (*_blossom_data)[left].next = arc;
1224 1243
      (*_blossom_data)[right].next = _graph.oppositeArc(arc);
1225 1244
    }
1226 1245

	
1227 1246
    void extendOnArc(const Arc& arc) {
1228 1247
      int base = _blossom_set->find(_graph.target(arc));
1229 1248
      int tree = _tree_set->find(base);
1230 1249

	
1231 1250
      int odd = _blossom_set->find(_graph.source(arc));
1232 1251
      _tree_set->insert(odd, tree);
1233 1252
      (*_blossom_data)[odd].status = ODD;
1234 1253
      matchedToOdd(odd);
1235 1254
      (*_blossom_data)[odd].pred = arc;
1236 1255

	
1237 1256
      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
1238 1257
      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1239 1258
      _tree_set->insert(even, tree);
1240 1259
      (*_blossom_data)[even].status = EVEN;
1241 1260
      matchedToEven(even, tree);
1242 1261
    }
1243 1262

	
1244 1263
    void shrinkOnEdge(const Edge& edge, int tree) {
1245 1264
      int nca = -1;
1246 1265
      std::vector<int> left_path, right_path;
1247 1266

	
1248 1267
      {
1249 1268
        std::set<int> left_set, right_set;
1250 1269
        int left = _blossom_set->find(_graph.u(edge));
1251 1270
        left_path.push_back(left);
1252 1271
        left_set.insert(left);
1253 1272

	
1254 1273
        int right = _blossom_set->find(_graph.v(edge));
1255 1274
        right_path.push_back(right);
1256 1275
        right_set.insert(right);
1257 1276

	
1258 1277
        while (true) {
1259 1278

	
1260 1279
          if ((*_blossom_data)[left].pred == INVALID) break;
1261 1280

	
1262 1281
          left =
1263 1282
            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1264 1283
          left_path.push_back(left);
1265 1284
          left =
1266 1285
            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1267 1286
          left_path.push_back(left);
1268 1287

	
1269 1288
          left_set.insert(left);
1270 1289

	
1271 1290
          if (right_set.find(left) != right_set.end()) {
1272 1291
            nca = left;
1273 1292
            break;
1274 1293
          }
1275 1294

	
1276 1295
          if ((*_blossom_data)[right].pred == INVALID) break;
1277 1296

	
1278 1297
          right =
1279 1298
            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1280 1299
          right_path.push_back(right);
1281 1300
          right =
1282 1301
            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1283 1302
          right_path.push_back(right);
1284 1303

	
1285 1304
          right_set.insert(right);
1286 1305

	
1287 1306
          if (left_set.find(right) != left_set.end()) {
1288 1307
            nca = right;
1289 1308
            break;
1290 1309
          }
1291 1310

	
1292 1311
        }
1293 1312

	
1294 1313
        if (nca == -1) {
1295 1314
          if ((*_blossom_data)[left].pred == INVALID) {
1296 1315
            nca = right;
1297 1316
            while (left_set.find(nca) == left_set.end()) {
1298 1317
              nca =
1299 1318
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1300 1319
              right_path.push_back(nca);
1301 1320
              nca =
1302 1321
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1303 1322
              right_path.push_back(nca);
1304 1323
            }
1305 1324
          } else {
1306 1325
            nca = left;
1307 1326
            while (right_set.find(nca) == right_set.end()) {
1308 1327
              nca =
1309 1328
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1310 1329
              left_path.push_back(nca);
1311 1330
              nca =
1312 1331
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1313 1332
              left_path.push_back(nca);
1314 1333
            }
1315 1334
          }
1316 1335
        }
1317 1336
      }
1318 1337

	
1319 1338
      std::vector<int> subblossoms;
1320 1339
      Arc prev;
1321 1340

	
1322 1341
      prev = _graph.direct(edge, true);
1323 1342
      for (int i = 0; left_path[i] != nca; i += 2) {
1324 1343
        subblossoms.push_back(left_path[i]);
1325 1344
        (*_blossom_data)[left_path[i]].next = prev;
1326 1345
        _tree_set->erase(left_path[i]);
1327 1346

	
1328 1347
        subblossoms.push_back(left_path[i + 1]);
1329 1348
        (*_blossom_data)[left_path[i + 1]].status = EVEN;
1330 1349
        oddToEven(left_path[i + 1], tree);
1331 1350
        _tree_set->erase(left_path[i + 1]);
1332 1351
        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
1333 1352
      }
1334 1353

	
1335 1354
      int k = 0;
1336 1355
      while (right_path[k] != nca) ++k;
1337 1356

	
1338 1357
      subblossoms.push_back(nca);
1339 1358
      (*_blossom_data)[nca].next = prev;
1340 1359

	
1341 1360
      for (int i = k - 2; i >= 0; i -= 2) {
1342 1361
        subblossoms.push_back(right_path[i + 1]);
1343 1362
        (*_blossom_data)[right_path[i + 1]].status = EVEN;
1344 1363
        oddToEven(right_path[i + 1], tree);
1345 1364
        _tree_set->erase(right_path[i + 1]);
1346 1365

	
1347 1366
        (*_blossom_data)[right_path[i + 1]].next =
1348 1367
          (*_blossom_data)[right_path[i + 1]].pred;
1349 1368

	
1350 1369
        subblossoms.push_back(right_path[i]);
1351 1370
        _tree_set->erase(right_path[i]);
1352 1371
      }
1353 1372

	
1354 1373
      int surface =
1355 1374
        _blossom_set->join(subblossoms.begin(), subblossoms.end());
1356 1375

	
1357 1376
      for (int i = 0; i < int(subblossoms.size()); ++i) {
1358 1377
        if (!_blossom_set->trivial(subblossoms[i])) {
1359 1378
          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1360 1379
        }
1361 1380
        (*_blossom_data)[subblossoms[i]].status = MATCHED;
1362 1381
      }
1363 1382

	
1364 1383
      (*_blossom_data)[surface].pot = -2 * _delta_sum;
1365 1384
      (*_blossom_data)[surface].offset = 0;
1366 1385
      (*_blossom_data)[surface].status = EVEN;
1367 1386
      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1368 1387
      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1369 1388

	
1370 1389
      _tree_set->insert(surface, tree);
1371 1390
      _tree_set->erase(nca);
1372 1391
    }
1373 1392

	
1374 1393
    void splitBlossom(int blossom) {
1375 1394
      Arc next = (*_blossom_data)[blossom].next;
1376 1395
      Arc pred = (*_blossom_data)[blossom].pred;
1377 1396

	
1378 1397
      int tree = _tree_set->find(blossom);
1379 1398

	
1380 1399
      (*_blossom_data)[blossom].status = MATCHED;
1381 1400
      oddToMatched(blossom);
1382 1401
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1383 1402
        _delta2->erase(blossom);
1384 1403
      }
1385 1404

	
1386 1405
      std::vector<int> subblossoms;
1387 1406
      _blossom_set->split(blossom, std::back_inserter(subblossoms));
1388 1407

	
1389 1408
      Value offset = (*_blossom_data)[blossom].offset;
1390 1409
      int b = _blossom_set->find(_graph.source(pred));
1391 1410
      int d = _blossom_set->find(_graph.source(next));
1392 1411

	
1393 1412
      int ib = -1, id = -1;
1394 1413
      for (int i = 0; i < int(subblossoms.size()); ++i) {
1395 1414
        if (subblossoms[i] == b) ib = i;
1396 1415
        if (subblossoms[i] == d) id = i;
1397 1416

	
1398 1417
        (*_blossom_data)[subblossoms[i]].offset = offset;
1399 1418
        if (!_blossom_set->trivial(subblossoms[i])) {
1400 1419
          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1401 1420
        }
1402 1421
        if (_blossom_set->classPrio(subblossoms[i]) !=
1403 1422
            std::numeric_limits<Value>::max()) {
1404 1423
          _delta2->push(subblossoms[i],
1405 1424
                        _blossom_set->classPrio(subblossoms[i]) -
1406 1425
                        (*_blossom_data)[subblossoms[i]].offset);
1407 1426
        }
1408 1427
      }
1409 1428

	
1410 1429
      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1411 1430
        for (int i = (id + 1) % subblossoms.size();
1412 1431
             i != ib; i = (i + 2) % subblossoms.size()) {
1413 1432
          int sb = subblossoms[i];
1414 1433
          int tb = subblossoms[(i + 1) % subblossoms.size()];
1415 1434
          (*_blossom_data)[sb].next =
1416 1435
            _graph.oppositeArc((*_blossom_data)[tb].next);
1417 1436
        }
1418 1437

	
1419 1438
        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1420 1439
          int sb = subblossoms[i];
1421 1440
          int tb = subblossoms[(i + 1) % subblossoms.size()];
1422 1441
          int ub = subblossoms[(i + 2) % subblossoms.size()];
1423 1442

	
1424 1443
          (*_blossom_data)[sb].status = ODD;
1425 1444
          matchedToOdd(sb);
1426 1445
          _tree_set->insert(sb, tree);
1427 1446
          (*_blossom_data)[sb].pred = pred;
1428 1447
          (*_blossom_data)[sb].next =
1429 1448
            _graph.oppositeArc((*_blossom_data)[tb].next);
1430 1449

	
1431 1450
          pred = (*_blossom_data)[ub].next;
1432 1451

	
1433 1452
          (*_blossom_data)[tb].status = EVEN;
1434 1453
          matchedToEven(tb, tree);
1435 1454
          _tree_set->insert(tb, tree);
1436 1455
          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1437 1456
        }
1438 1457

	
1439 1458
        (*_blossom_data)[subblossoms[id]].status = ODD;
1440 1459
        matchedToOdd(subblossoms[id]);
1441 1460
        _tree_set->insert(subblossoms[id], tree);
1442 1461
        (*_blossom_data)[subblossoms[id]].next = next;
1443 1462
        (*_blossom_data)[subblossoms[id]].pred = pred;
1444 1463

	
1445 1464
      } else {
1446 1465

	
1447 1466
        for (int i = (ib + 1) % subblossoms.size();
1448 1467
             i != id; i = (i + 2) % subblossoms.size()) {
1449 1468
          int sb = subblossoms[i];
1450 1469
          int tb = subblossoms[(i + 1) % subblossoms.size()];
1451 1470
          (*_blossom_data)[sb].next =
1452 1471
            _graph.oppositeArc((*_blossom_data)[tb].next);
1453 1472
        }
1454 1473

	
1455 1474
        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1456 1475
          int sb = subblossoms[i];
1457 1476
          int tb = subblossoms[(i + 1) % subblossoms.size()];
1458 1477
          int ub = subblossoms[(i + 2) % subblossoms.size()];
1459 1478

	
1460 1479
          (*_blossom_data)[sb].status = ODD;
1461 1480
          matchedToOdd(sb);
1462 1481
          _tree_set->insert(sb, tree);
1463 1482
          (*_blossom_data)[sb].next = next;
1464 1483
          (*_blossom_data)[sb].pred =
1465 1484
            _graph.oppositeArc((*_blossom_data)[tb].next);
1466 1485

	
1467 1486
          (*_blossom_data)[tb].status = EVEN;
1468 1487
          matchedToEven(tb, tree);
1469 1488
          _tree_set->insert(tb, tree);
1470 1489
          (*_blossom_data)[tb].pred =
1471 1490
            (*_blossom_data)[tb].next =
1472 1491
            _graph.oppositeArc((*_blossom_data)[ub].next);
1473 1492
          next = (*_blossom_data)[ub].next;
1474 1493
        }
1475 1494

	
1476 1495
        (*_blossom_data)[subblossoms[ib]].status = ODD;
1477 1496
        matchedToOdd(subblossoms[ib]);
1478 1497
        _tree_set->insert(subblossoms[ib], tree);
1479 1498
        (*_blossom_data)[subblossoms[ib]].next = next;
1480 1499
        (*_blossom_data)[subblossoms[ib]].pred = pred;
1481 1500
      }
1482 1501
      _tree_set->erase(blossom);
1483 1502
    }
1484 1503

	
1485 1504
    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1486 1505
      if (_blossom_set->trivial(blossom)) {
1487 1506
        int bi = (*_node_index)[base];
1488 1507
        Value pot = (*_node_data)[bi].pot;
1489 1508

	
1490 1509
        (*_matching)[base] = matching;
1491 1510
        _blossom_node_list.push_back(base);
1492 1511
        (*_node_potential)[base] = pot;
1493 1512
      } else {
1494 1513

	
1495 1514
        Value pot = (*_blossom_data)[blossom].pot;
1496 1515
        int bn = _blossom_node_list.size();
1497 1516

	
1498 1517
        std::vector<int> subblossoms;
1499 1518
        _blossom_set->split(blossom, std::back_inserter(subblossoms));
1500 1519
        int b = _blossom_set->find(base);
1501 1520
        int ib = -1;
1502 1521
        for (int i = 0; i < int(subblossoms.size()); ++i) {
1503 1522
          if (subblossoms[i] == b) { ib = i; break; }
1504 1523
        }
1505 1524

	
1506 1525
        for (int i = 1; i < int(subblossoms.size()); i += 2) {
1507 1526
          int sb = subblossoms[(ib + i) % subblossoms.size()];
1508 1527
          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1509 1528

	
1510 1529
          Arc m = (*_blossom_data)[tb].next;
1511 1530
          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1512 1531
          extractBlossom(tb, _graph.source(m), m);
1513 1532
        }
1514 1533
        extractBlossom(subblossoms[ib], base, matching);
1515 1534

	
1516 1535
        int en = _blossom_node_list.size();
1517 1536

	
1518 1537
        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1519 1538
      }
1520 1539
    }
1521 1540

	
1522 1541
    void extractMatching() {
1523 1542
      std::vector<int> blossoms;
1524 1543
      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1525 1544
        blossoms.push_back(c);
1526 1545
      }
1527 1546

	
1528 1547
      for (int i = 0; i < int(blossoms.size()); ++i) {
1529 1548
        if ((*_blossom_data)[blossoms[i]].next != INVALID) {
1530 1549

	
1531 1550
          Value offset = (*_blossom_data)[blossoms[i]].offset;
1532 1551
          (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1533 1552
          for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1534 1553
               n != INVALID; ++n) {
1535 1554
            (*_node_data)[(*_node_index)[n]].pot -= offset;
1536 1555
          }
1537 1556

	
1538 1557
          Arc matching = (*_blossom_data)[blossoms[i]].next;
1539 1558
          Node base = _graph.source(matching);
1540 1559
          extractBlossom(blossoms[i], base, matching);
1541 1560
        } else {
1542 1561
          Node base = (*_blossom_data)[blossoms[i]].base;
1543 1562
          extractBlossom(blossoms[i], base, INVALID);
1544 1563
        }
1545 1564
      }
1546 1565
    }
1547 1566

	
1548 1567
  public:
1549 1568

	
1550 1569
    /// \brief Constructor
1551 1570
    ///
1552 1571
    /// Constructor.
1553 1572
    MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1554 1573
      : _graph(graph), _weight(weight), _matching(0),
1555 1574
        _node_potential(0), _blossom_potential(), _blossom_node_list(),
1556 1575
        _node_num(0), _blossom_num(0),
1557 1576

	
1558 1577
        _blossom_index(0), _blossom_set(0), _blossom_data(0),
1559 1578
        _node_index(0), _node_heap_index(0), _node_data(0),
1560 1579
        _tree_set_index(0), _tree_set(0),
1561 1580

	
1562 1581
        _delta1_index(0), _delta1(0),
1563 1582
        _delta2_index(0), _delta2(0),
1564 1583
        _delta3_index(0), _delta3(0),
1565 1584
        _delta4_index(0), _delta4(0),
1566 1585

	
1567 1586
        _delta_sum(), _unmatched(0),
1568 1587

	
1569 1588
        _fractional(0)
1570 1589
    {}
1571 1590

	
1572 1591
    ~MaxWeightedMatching() {
1573 1592
      destroyStructures();
1574 1593
      if (_fractional) {
1575 1594
        delete _fractional;
1576 1595
      }
1577 1596
    }
1578 1597

	
1579 1598
    /// \name Execution Control
1580 1599
    /// The simplest way to execute the algorithm is to use the
1581 1600
    /// \ref run() member function.
1582 1601

	
1583 1602
    ///@{
1584 1603

	
1585 1604
    /// \brief Initialize the algorithm
1586 1605
    ///
1587 1606
    /// This function initializes the algorithm.
1588 1607
    void init() {
1589 1608
      createStructures();
1590 1609

	
1610
      _blossom_node_list.clear();
1611
      _blossom_potential.clear();
1612

	
1591 1613
      for (ArcIt e(_graph); e != INVALID; ++e) {
1592 1614
        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1593 1615
      }
1594 1616
      for (NodeIt n(_graph); n != INVALID; ++n) {
1595 1617
        (*_delta1_index)[n] = _delta1->PRE_HEAP;
1596 1618
      }
1597 1619
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1598 1620
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1599 1621
      }
1600 1622
      for (int i = 0; i < _blossom_num; ++i) {
1601 1623
        (*_delta2_index)[i] = _delta2->PRE_HEAP;
1602 1624
        (*_delta4_index)[i] = _delta4->PRE_HEAP;
1603 1625
      }
1604

	
1626
      
1605 1627
      _unmatched = _node_num;
1606 1628

	
1629
      _delta1->clear();
1630
      _delta2->clear();
1631
      _delta3->clear();
1632
      _delta4->clear();
1633
      _blossom_set->clear();
1634
      _tree_set->clear();
1635

	
1607 1636
      int index = 0;
1608 1637
      for (NodeIt n(_graph); n != INVALID; ++n) {
1609 1638
        Value max = 0;
1610 1639
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1611 1640
          if (_graph.target(e) == n) continue;
1612 1641
          if ((dualScale * _weight[e]) / 2 > max) {
1613 1642
            max = (dualScale * _weight[e]) / 2;
1614 1643
          }
1615 1644
        }
1616 1645
        (*_node_index)[n] = index;
1646
        (*_node_data)[index].heap_index.clear();
1647
        (*_node_data)[index].heap.clear();
1617 1648
        (*_node_data)[index].pot = max;
1618 1649
        _delta1->push(n, max);
1619 1650
        int blossom =
1620 1651
          _blossom_set->insert(n, std::numeric_limits<Value>::max());
1621 1652

	
1622 1653
        _tree_set->insert(blossom);
1623 1654

	
1624 1655
        (*_blossom_data)[blossom].status = EVEN;
1625 1656
        (*_blossom_data)[blossom].pred = INVALID;
1626 1657
        (*_blossom_data)[blossom].next = INVALID;
1627 1658
        (*_blossom_data)[blossom].pot = 0;
1628 1659
        (*_blossom_data)[blossom].offset = 0;
1629 1660
        ++index;
1630 1661
      }
1631 1662
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1632 1663
        int si = (*_node_index)[_graph.u(e)];
1633 1664
        int ti = (*_node_index)[_graph.v(e)];
1634 1665
        if (_graph.u(e) != _graph.v(e)) {
1635 1666
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1636 1667
                            dualScale * _weight[e]) / 2);
1637 1668
        }
1638 1669
      }
1639 1670
    }
1640 1671

	
1641 1672
    /// \brief Initialize the algorithm with fractional matching
1642 1673
    ///
1643 1674
    /// This function initializes the algorithm with a fractional
1644 1675
    /// matching. This initialization is also called jumpstart heuristic.
1645 1676
    void fractionalInit() {
1646 1677
      createStructures();
1647 1678
      
1648 1679
      if (_fractional == 0) {
1649 1680
        _fractional = new FractionalMatching(_graph, _weight, false);
1650 1681
      }
1651 1682
      _fractional->run();
1652 1683

	
1653 1684
      for (ArcIt e(_graph); e != INVALID; ++e) {
1654 1685
        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1655 1686
      }
1656 1687
      for (NodeIt n(_graph); n != INVALID; ++n) {
1657 1688
        (*_delta1_index)[n] = _delta1->PRE_HEAP;
1658 1689
      }
1659 1690
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1660 1691
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1661 1692
      }
1662 1693
      for (int i = 0; i < _blossom_num; ++i) {
1663 1694
        (*_delta2_index)[i] = _delta2->PRE_HEAP;
1664 1695
        (*_delta4_index)[i] = _delta4->PRE_HEAP;
1665 1696
      }
1666 1697

	
1667 1698
      _unmatched = 0;
1668 1699

	
1669 1700
      int index = 0;
1670 1701
      for (NodeIt n(_graph); n != INVALID; ++n) {
1671 1702
        Value pot = _fractional->nodeValue(n);
1672 1703
        (*_node_index)[n] = index;
1673 1704
        (*_node_data)[index].pot = pot;
1674 1705
        int blossom =
1675 1706
          _blossom_set->insert(n, std::numeric_limits<Value>::max());
1676 1707

	
1677 1708
        (*_blossom_data)[blossom].status = MATCHED;
1678 1709
        (*_blossom_data)[blossom].pred = INVALID;
1679 1710
        (*_blossom_data)[blossom].next = _fractional->matching(n);
1680 1711
        if (_fractional->matching(n) == INVALID) {
1681 1712
          (*_blossom_data)[blossom].base = n;
1682 1713
        }
1683 1714
        (*_blossom_data)[blossom].pot = 0;
1684 1715
        (*_blossom_data)[blossom].offset = 0;
1685 1716
        ++index;
1686 1717
      }
1687 1718

	
1688 1719
      typename Graph::template NodeMap<bool> processed(_graph, false);
1689 1720
      for (NodeIt n(_graph); n != INVALID; ++n) {
1690 1721
        if (processed[n]) continue;
1691 1722
        processed[n] = true;
1692 1723
        if (_fractional->matching(n) == INVALID) continue;
1693 1724
        int num = 1;
1694 1725
        Node v = _graph.target(_fractional->matching(n));
1695 1726
        while (n != v) {
1696 1727
          processed[v] = true;
1697 1728
          v = _graph.target(_fractional->matching(v));
1698 1729
          ++num;
1699 1730
        }
1700 1731

	
1701 1732
        if (num % 2 == 1) {
1702 1733
          std::vector<int> subblossoms(num);
1703 1734

	
1704 1735
          subblossoms[--num] = _blossom_set->find(n);
1705 1736
          _delta1->push(n, _fractional->nodeValue(n));
1706 1737
          v = _graph.target(_fractional->matching(n));
1707 1738
          while (n != v) {
1708 1739
            subblossoms[--num] = _blossom_set->find(v);
1709 1740
            _delta1->push(v, _fractional->nodeValue(v));
1710 1741
            v = _graph.target(_fractional->matching(v));            
1711 1742
          }
1712 1743
          
1713 1744
          int surface = 
1714 1745
            _blossom_set->join(subblossoms.begin(), subblossoms.end());
1715 1746
          (*_blossom_data)[surface].status = EVEN;
1716 1747
          (*_blossom_data)[surface].pred = INVALID;
1717 1748
          (*_blossom_data)[surface].next = INVALID;
1718 1749
          (*_blossom_data)[surface].pot = 0;
1719 1750
          (*_blossom_data)[surface].offset = 0;
1720 1751
          
1721 1752
          _tree_set->insert(surface);
1722 1753
          ++_unmatched;
1723 1754
        }
1724 1755
      }
1725 1756

	
1726 1757
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1727 1758
        int si = (*_node_index)[_graph.u(e)];
1728 1759
        int sb = _blossom_set->find(_graph.u(e));
1729 1760
        int ti = (*_node_index)[_graph.v(e)];
1730 1761
        int tb = _blossom_set->find(_graph.v(e));
1731 1762
        if ((*_blossom_data)[sb].status == EVEN &&
1732 1763
            (*_blossom_data)[tb].status == EVEN && sb != tb) {
1733 1764
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1734 1765
                            dualScale * _weight[e]) / 2);
1735 1766
        }
1736 1767
      }
1737 1768

	
1738 1769
      for (NodeIt n(_graph); n != INVALID; ++n) {
1739 1770
        int nb = _blossom_set->find(n);
1740 1771
        if ((*_blossom_data)[nb].status != MATCHED) continue;
1741 1772
        int ni = (*_node_index)[n];
1742 1773

	
1743 1774
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1744 1775
          Node v = _graph.target(e);
1745 1776
          int vb = _blossom_set->find(v);
1746 1777
          int vi = (*_node_index)[v];
1747 1778

	
1748 1779
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1749 1780
            dualScale * _weight[e];
1750 1781

	
1751 1782
          if ((*_blossom_data)[vb].status == EVEN) {
1752 1783

	
1753 1784
            int vt = _tree_set->find(vb);
1754 1785

	
1755 1786
            typename std::map<int, Arc>::iterator it =
1756 1787
              (*_node_data)[ni].heap_index.find(vt);
1757 1788

	
1758 1789
            if (it != (*_node_data)[ni].heap_index.end()) {
1759 1790
              if ((*_node_data)[ni].heap[it->second] > rw) {
1760 1791
                (*_node_data)[ni].heap.replace(it->second, e);
1761 1792
                (*_node_data)[ni].heap.decrease(e, rw);
1762 1793
                it->second = e;
1763 1794
              }
1764 1795
            } else {
1765 1796
              (*_node_data)[ni].heap.push(e, rw);
1766 1797
              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
1767 1798
            }
1768 1799
          }
1769 1800
        }
1770 1801
            
1771 1802
        if (!(*_node_data)[ni].heap.empty()) {
1772 1803
          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1773 1804
          _delta2->push(nb, _blossom_set->classPrio(nb));
1774 1805
        }
1775 1806
      }
1776 1807
    }
1777 1808

	
1778 1809
    /// \brief Start the algorithm
1779 1810
    ///
1780 1811
    /// This function starts the algorithm.
1781 1812
    ///
1782 1813
    /// \pre \ref init() or \ref fractionalInit() must be called
1783 1814
    /// before using this function.
1784 1815
    void start() {
1785 1816
      enum OpType {
1786 1817
        D1, D2, D3, D4
1787 1818
      };
1788 1819

	
1789 1820
      while (_unmatched > 0) {
1790 1821
        Value d1 = !_delta1->empty() ?
1791 1822
          _delta1->prio() : std::numeric_limits<Value>::max();
1792 1823

	
1793 1824
        Value d2 = !_delta2->empty() ?
1794 1825
          _delta2->prio() : std::numeric_limits<Value>::max();
1795 1826

	
1796 1827
        Value d3 = !_delta3->empty() ?
1797 1828
          _delta3->prio() : std::numeric_limits<Value>::max();
1798 1829

	
1799 1830
        Value d4 = !_delta4->empty() ?
1800 1831
          _delta4->prio() : std::numeric_limits<Value>::max();
1801 1832

	
1802 1833
        _delta_sum = d3; OpType ot = D3;
1803 1834
        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1804 1835
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1805 1836
        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1806 1837

	
1807 1838
        switch (ot) {
1808 1839
        case D1:
1809 1840
          {
1810 1841
            Node n = _delta1->top();
1811 1842
            unmatchNode(n);
1812 1843
            --_unmatched;
1813 1844
          }
1814 1845
          break;
1815 1846
        case D2:
1816 1847
          {
1817 1848
            int blossom = _delta2->top();
1818 1849
            Node n = _blossom_set->classTop(blossom);
1819 1850
            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
1820 1851
            if ((*_blossom_data)[blossom].next == INVALID) {
1821 1852
              augmentOnArc(a);
1822 1853
              --_unmatched;
1823 1854
            } else {
1824 1855
              extendOnArc(a);
1825 1856
            }
1826 1857
          }
1827 1858
          break;
1828 1859
        case D3:
1829 1860
          {
1830 1861
            Edge e = _delta3->top();
1831 1862

	
1832 1863
            int left_blossom = _blossom_set->find(_graph.u(e));
1833 1864
            int right_blossom = _blossom_set->find(_graph.v(e));
1834 1865

	
1835 1866
            if (left_blossom == right_blossom) {
1836 1867
              _delta3->pop();
1837 1868
            } else {
1838 1869
              int left_tree = _tree_set->find(left_blossom);
1839 1870
              int right_tree = _tree_set->find(right_blossom);
1840 1871

	
1841 1872
              if (left_tree == right_tree) {
1842 1873
                shrinkOnEdge(e, left_tree);
1843 1874
              } else {
1844 1875
                augmentOnEdge(e);
1845 1876
                _unmatched -= 2;
1846 1877
              }
1847 1878
            }
1848 1879
          } break;
1849 1880
        case D4:
1850 1881
          splitBlossom(_delta4->top());
1851 1882
          break;
1852 1883
        }
1853 1884
      }
1854 1885
      extractMatching();
1855 1886
    }
1856 1887

	
1857 1888
    /// \brief Run the algorithm.
1858 1889
    ///
1859 1890
    /// This method runs the \c %MaxWeightedMatching algorithm.
1860 1891
    ///
1861 1892
    /// \note mwm.run() is just a shortcut of the following code.
1862 1893
    /// \code
1863 1894
    ///   mwm.fractionalInit();
1864 1895
    ///   mwm.start();
1865 1896
    /// \endcode
1866 1897
    void run() {
1867 1898
      fractionalInit();
1868 1899
      start();
1869 1900
    }
1870 1901

	
1871 1902
    /// @}
1872 1903

	
1873 1904
    /// \name Primal Solution
1874 1905
    /// Functions to get the primal solution, i.e. the maximum weighted
1875 1906
    /// matching.\n
1876 1907
    /// Either \ref run() or \ref start() function should be called before
1877 1908
    /// using them.
1878 1909

	
1879 1910
    /// @{
1880 1911

	
1881 1912
    /// \brief Return the weight of the matching.
1882 1913
    ///
1883 1914
    /// This function returns the weight of the found matching.
1884 1915
    ///
1885 1916
    /// \pre Either run() or start() must be called before using this function.
1886 1917
    Value matchingWeight() const {
1887 1918
      Value sum = 0;
1888 1919
      for (NodeIt n(_graph); n != INVALID; ++n) {
1889 1920
        if ((*_matching)[n] != INVALID) {
1890 1921
          sum += _weight[(*_matching)[n]];
1891 1922
        }
1892 1923
      }
1893 1924
      return sum / 2;
1894 1925
    }
1895 1926

	
1896 1927
    /// \brief Return the size (cardinality) of the matching.
1897 1928
    ///
1898 1929
    /// This function returns the size (cardinality) of the found matching.
1899 1930
    ///
1900 1931
    /// \pre Either run() or start() must be called before using this function.
1901 1932
    int matchingSize() const {
1902 1933
      int num = 0;
1903 1934
      for (NodeIt n(_graph); n != INVALID; ++n) {
1904 1935
        if ((*_matching)[n] != INVALID) {
1905 1936
          ++num;
1906 1937
        }
1907 1938
      }
1908 1939
      return num /= 2;
1909 1940
    }
1910 1941

	
1911 1942
    /// \brief Return \c true if the given edge is in the matching.
1912 1943
    ///
1913 1944
    /// This function returns \c true if the given edge is in the found
1914 1945
    /// matching.
1915 1946
    ///
1916 1947
    /// \pre Either run() or start() must be called before using this function.
1917 1948
    bool matching(const Edge& edge) const {
1918 1949
      return edge == (*_matching)[_graph.u(edge)];
1919 1950
    }
1920 1951

	
1921 1952
    /// \brief Return the matching arc (or edge) incident to the given node.
1922 1953
    ///
1923 1954
    /// This function returns the matching arc (or edge) incident to the
1924 1955
    /// given node in the found matching or \c INVALID if the node is
1925 1956
    /// not covered by the matching.
1926 1957
    ///
1927 1958
    /// \pre Either run() or start() must be called before using this function.
1928 1959
    Arc matching(const Node& node) const {
1929 1960
      return (*_matching)[node];
1930 1961
    }
1931 1962

	
1932 1963
    /// \brief Return a const reference to the matching map.
1933 1964
    ///
1934 1965
    /// This function returns a const reference to a node map that stores
1935 1966
    /// the matching arc (or edge) incident to each node.
1936 1967
    const MatchingMap& matchingMap() const {
1937 1968
      return *_matching;
1938 1969
    }
1939 1970

	
1940 1971
    /// \brief Return the mate of the given node.
1941 1972
    ///
1942 1973
    /// This function returns the mate of the given node in the found
1943 1974
    /// matching or \c INVALID if the node is not covered by the matching.
1944 1975
    ///
1945 1976
    /// \pre Either run() or start() must be called before using this function.
1946 1977
    Node mate(const Node& node) const {
1947 1978
      return (*_matching)[node] != INVALID ?
1948 1979
        _graph.target((*_matching)[node]) : INVALID;
1949 1980
    }
1950 1981

	
1951 1982
    /// @}
1952 1983

	
1953 1984
    /// \name Dual Solution
1954 1985
    /// Functions to get the dual solution.\n
1955 1986
    /// Either \ref run() or \ref start() function should be called before
1956 1987
    /// using them.
1957 1988

	
1958 1989
    /// @{
1959 1990

	
1960 1991
    /// \brief Return the value of the dual solution.
1961 1992
    ///
1962 1993
    /// This function returns the value of the dual solution.
1963 1994
    /// It should be equal to the primal value scaled by \ref dualScale
1964 1995
    /// "dual scale".
1965 1996
    ///
1966 1997
    /// \pre Either run() or start() must be called before using this function.
1967 1998
    Value dualValue() const {
1968 1999
      Value sum = 0;
1969 2000
      for (NodeIt n(_graph); n != INVALID; ++n) {
1970 2001
        sum += nodeValue(n);
1971 2002
      }
1972 2003
      for (int i = 0; i < blossomNum(); ++i) {
1973 2004
        sum += blossomValue(i) * (blossomSize(i) / 2);
1974 2005
      }
1975 2006
      return sum;
1976 2007
    }
1977 2008

	
1978 2009
    /// \brief Return the dual value (potential) of the given node.
1979 2010
    ///
1980 2011
    /// This function returns the dual value (potential) of the given node.
1981 2012
    ///
1982 2013
    /// \pre Either run() or start() must be called before using this function.
1983 2014
    Value nodeValue(const Node& n) const {
1984 2015
      return (*_node_potential)[n];
1985 2016
    }
1986 2017

	
1987 2018
    /// \brief Return the number of the blossoms in the basis.
1988 2019
    ///
1989 2020
    /// This function returns the number of the blossoms in the basis.
1990 2021
    ///
1991 2022
    /// \pre Either run() or start() must be called before using this function.
1992 2023
    /// \see BlossomIt
1993 2024
    int blossomNum() const {
1994 2025
      return _blossom_potential.size();
1995 2026
    }
1996 2027

	
1997 2028
    /// \brief Return the number of the nodes in the given blossom.
1998 2029
    ///
1999 2030
    /// This function returns the number of the nodes in the given blossom.
2000 2031
    ///
2001 2032
    /// \pre Either run() or start() must be called before using this function.
2002 2033
    /// \see BlossomIt
2003 2034
    int blossomSize(int k) const {
2004 2035
      return _blossom_potential[k].end - _blossom_potential[k].begin;
2005 2036
    }
2006 2037

	
2007 2038
    /// \brief Return the dual value (ptential) of the given blossom.
2008 2039
    ///
2009 2040
    /// This function returns the dual value (ptential) of the given blossom.
2010 2041
    ///
2011 2042
    /// \pre Either run() or start() must be called before using this function.
2012 2043
    Value blossomValue(int k) const {
2013 2044
      return _blossom_potential[k].value;
2014 2045
    }
2015 2046

	
2016 2047
    /// \brief Iterator for obtaining the nodes of a blossom.
2017 2048
    ///
2018 2049
    /// This class provides an iterator for obtaining the nodes of the
2019 2050
    /// given blossom. It lists a subset of the nodes.
2020 2051
    /// Before using this iterator, you must allocate a
2021 2052
    /// MaxWeightedMatching class and execute it.
2022 2053
    class BlossomIt {
2023 2054
    public:
2024 2055

	
2025 2056
      /// \brief Constructor.
2026 2057
      ///
2027 2058
      /// Constructor to get the nodes of the given variable.
2028 2059
      ///
2029 2060
      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
2030 2061
      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
2031 2062
      /// called before initializing this iterator.
2032 2063
      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
2033 2064
        : _algorithm(&algorithm)
2034 2065
      {
2035 2066
        _index = _algorithm->_blossom_potential[variable].begin;
2036 2067
        _last = _algorithm->_blossom_potential[variable].end;
2037 2068
      }
2038 2069

	
2039 2070
      /// \brief Conversion to \c Node.
2040 2071
      ///
2041 2072
      /// Conversion to \c Node.
2042 2073
      operator Node() const {
2043 2074
        return _algorithm->_blossom_node_list[_index];
2044 2075
      }
2045 2076

	
2046 2077
      /// \brief Increment operator.
2047 2078
      ///
2048 2079
      /// Increment operator.
2049 2080
      BlossomIt& operator++() {
2050 2081
        ++_index;
2051 2082
        return *this;
2052 2083
      }
2053 2084

	
2054 2085
      /// \brief Validity checking
2055 2086
      ///
2056 2087
      /// Checks whether the iterator is invalid.
2057 2088
      bool operator==(Invalid) const { return _index == _last; }
2058 2089

	
2059 2090
      /// \brief Validity checking
2060 2091
      ///
2061 2092
      /// Checks whether the iterator is valid.
2062 2093
      bool operator!=(Invalid) const { return _index != _last; }
2063 2094

	
2064 2095
    private:
2065 2096
      const MaxWeightedMatching* _algorithm;
2066 2097
      int _last;
2067 2098
      int _index;
2068 2099
    };
2069 2100

	
2070 2101
    /// @}
2071 2102

	
2072 2103
  };
2073 2104

	
2074 2105
  /// \ingroup matching
2075 2106
  ///
2076 2107
  /// \brief Weighted perfect matching in general graphs
2077 2108
  ///
2078 2109
  /// This class provides an efficient implementation of Edmond's
2079 2110
  /// maximum weighted perfect matching algorithm. The implementation
2080 2111
  /// is based on extensive use of priority queues and provides
2081 2112
  /// \f$O(nm\log n)\f$ time complexity.
2082 2113
  ///
2083 2114
  /// The maximum weighted perfect matching problem is to find a subset of
2084 2115
  /// the edges in an undirected graph with maximum overall weight for which
2085 2116
  /// each node has exactly one incident edge.
2086 2117
  /// It can be formulated with the following linear program.
2087 2118
  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
2088 2119
  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
2089 2120
      \quad \forall B\in\mathcal{O}\f] */
2090 2121
  /// \f[x_e \ge 0\quad \forall e\in E\f]
2091 2122
  /// \f[\max \sum_{e\in E}x_ew_e\f]
2092 2123
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
2093 2124
  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
2094 2125
  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
2095 2126
  /// subsets of the nodes.
2096 2127
  ///
2097 2128
  /// The algorithm calculates an optimal matching and a proof of the
2098 2129
  /// optimality. The solution of the dual problem can be used to check
2099 2130
  /// the result of the algorithm. The dual linear problem is the
2100 2131
  /// following.
2101 2132
  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
2102 2133
      w_{uv} \quad \forall uv\in E\f] */
2103 2134
  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
2104 2135
  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
2105 2136
      \frac{\vert B \vert - 1}{2}z_B\f] */
2106 2137
  ///
2107 2138
  /// The algorithm can be executed with the run() function.
2108 2139
  /// After it the matching (the primal solution) and the dual solution
2109 2140
  /// can be obtained using the query functions and the
2110 2141
  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
2111 2142
  /// which is able to iterate on the nodes of a blossom.
2112 2143
  /// If the value type is integer, then the dual solution is multiplied
2113 2144
  /// by \ref MaxWeightedMatching::dualScale "4".
2114 2145
  ///
2115 2146
  /// \tparam GR The undirected graph type the algorithm runs on.
2116 2147
  /// \tparam WM The type edge weight map. The default type is
2117 2148
  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
2118 2149
#ifdef DOXYGEN
2119 2150
  template <typename GR, typename WM>
2120 2151
#else
2121 2152
  template <typename GR,
2122 2153
            typename WM = typename GR::template EdgeMap<int> >
2123 2154
#endif
2124 2155
  class MaxWeightedPerfectMatching {
2125 2156
  public:
2126 2157

	
2127 2158
    /// The graph type of the algorithm
2128 2159
    typedef GR Graph;
2129 2160
    /// The type of the edge weight map
2130 2161
    typedef WM WeightMap;
2131 2162
    /// The value type of the edge weights
2132 2163
    typedef typename WeightMap::Value Value;
2133 2164

	
2134 2165
    /// \brief Scaling factor for dual solution
2135 2166
    ///
2136 2167
    /// Scaling factor for dual solution, it is equal to 4 or 1
2137 2168
    /// according to the value type.
2138 2169
    static const int dualScale =
2139 2170
      std::numeric_limits<Value>::is_integer ? 4 : 1;
2140 2171

	
2141 2172
    /// The type of the matching map
2142 2173
    typedef typename Graph::template NodeMap<typename Graph::Arc>
2143 2174
    MatchingMap;
2144 2175

	
2145 2176
  private:
2146 2177

	
2147 2178
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
2148 2179

	
2149 2180
    typedef typename Graph::template NodeMap<Value> NodePotential;
2150 2181
    typedef std::vector<Node> BlossomNodeList;
2151 2182

	
2152 2183
    struct BlossomVariable {
2153 2184
      int begin, end;
2154 2185
      Value value;
2155 2186

	
2156 2187
      BlossomVariable(int _begin, int _end, Value _value)
2157 2188
        : begin(_begin), end(_end), value(_value) {}
2158 2189

	
2159 2190
    };
2160 2191

	
2161 2192
    typedef std::vector<BlossomVariable> BlossomPotential;
2162 2193

	
2163 2194
    const Graph& _graph;
2164 2195
    const WeightMap& _weight;
2165 2196

	
2166 2197
    MatchingMap* _matching;
2167 2198

	
2168 2199
    NodePotential* _node_potential;
2169 2200

	
2170 2201
    BlossomPotential _blossom_potential;
2171 2202
    BlossomNodeList _blossom_node_list;
2172 2203

	
2173 2204
    int _node_num;
2174 2205
    int _blossom_num;
2175 2206

	
2176 2207
    typedef RangeMap<int> IntIntMap;
2177 2208

	
2178 2209
    enum Status {
2179 2210
      EVEN = -1, MATCHED = 0, ODD = 1
2180 2211
    };
2181 2212

	
2182 2213
    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
2183 2214
    struct BlossomData {
2184 2215
      int tree;
2185 2216
      Status status;
2186 2217
      Arc pred, next;
2187 2218
      Value pot, offset;
2188 2219
    };
2189 2220

	
2190 2221
    IntNodeMap *_blossom_index;
2191 2222
    BlossomSet *_blossom_set;
2192 2223
    RangeMap<BlossomData>* _blossom_data;
2193 2224

	
2194 2225
    IntNodeMap *_node_index;
2195 2226
    IntArcMap *_node_heap_index;
2196 2227

	
2197 2228
    struct NodeData {
2198 2229

	
2199 2230
      NodeData(IntArcMap& node_heap_index)
2200 2231
        : heap(node_heap_index) {}
2201 2232

	
2202 2233
      int blossom;
2203 2234
      Value pot;
2204 2235
      BinHeap<Value, IntArcMap> heap;
2205 2236
      std::map<int, Arc> heap_index;
2206 2237

	
2207 2238
      int tree;
2208 2239
    };
2209 2240

	
2210 2241
    RangeMap<NodeData>* _node_data;
2211 2242

	
2212 2243
    typedef ExtendFindEnum<IntIntMap> TreeSet;
2213 2244

	
2214 2245
    IntIntMap *_tree_set_index;
2215 2246
    TreeSet *_tree_set;
2216 2247

	
2217 2248
    IntIntMap *_delta2_index;
2218 2249
    BinHeap<Value, IntIntMap> *_delta2;
2219 2250

	
2220 2251
    IntEdgeMap *_delta3_index;
2221 2252
    BinHeap<Value, IntEdgeMap> *_delta3;
2222 2253

	
2223 2254
    IntIntMap *_delta4_index;
2224 2255
    BinHeap<Value, IntIntMap> *_delta4;
2225 2256

	
2226 2257
    Value _delta_sum;
2227 2258
    int _unmatched;
2228 2259

	
2229 2260
    typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap> 
2230 2261
    FractionalMatching;
2231 2262
    FractionalMatching *_fractional;
2232 2263

	
2233 2264
    void createStructures() {
2234 2265
      _node_num = countNodes(_graph);
2235 2266
      _blossom_num = _node_num * 3 / 2;
2236 2267

	
2237 2268
      if (!_matching) {
2238 2269
        _matching = new MatchingMap(_graph);
2239 2270
      }
2271

	
2240 2272
      if (!_node_potential) {
2241 2273
        _node_potential = new NodePotential(_graph);
2242 2274
      }
2275

	
2243 2276
      if (!_blossom_set) {
2244 2277
        _blossom_index = new IntNodeMap(_graph);
2245 2278
        _blossom_set = new BlossomSet(*_blossom_index);
2246 2279
        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2280
      } else if (_blossom_data->size() != _blossom_num) {
2281
        delete _blossom_data;
2282
        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2247 2283
      }
2248 2284

	
2249 2285
      if (!_node_index) {
2250 2286
        _node_index = new IntNodeMap(_graph);
2251 2287
        _node_heap_index = new IntArcMap(_graph);
2252 2288
        _node_data = new RangeMap<NodeData>(_node_num,
2253 2289
                                            NodeData(*_node_heap_index));
2290
      } else if (_node_data->size() != _node_num) {
2291
        delete _node_data;
2292
        _node_data = new RangeMap<NodeData>(_node_num,
2293
                                            NodeData(*_node_heap_index));
2254 2294
      }
2255 2295

	
2256 2296
      if (!_tree_set) {
2257 2297
        _tree_set_index = new IntIntMap(_blossom_num);
2258 2298
        _tree_set = new TreeSet(*_tree_set_index);
2299
      } else {
2300
        _tree_set_index->resize(_blossom_num);
2259 2301
      }
2302

	
2260 2303
      if (!_delta2) {
2261 2304
        _delta2_index = new IntIntMap(_blossom_num);
2262 2305
        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2306
      } else {
2307
        _delta2_index->resize(_blossom_num);
2263 2308
      }
2309

	
2264 2310
      if (!_delta3) {
2265 2311
        _delta3_index = new IntEdgeMap(_graph);
2266 2312
        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
2267 2313
      }
2314

	
2268 2315
      if (!_delta4) {
2269 2316
        _delta4_index = new IntIntMap(_blossom_num);
2270 2317
        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2318
      } else {
2319
        _delta4_index->resize(_blossom_num);
2271 2320
      }
2272 2321
    }
2273 2322

	
2274 2323
    void destroyStructures() {
2275 2324
      if (_matching) {
2276 2325
        delete _matching;
2277 2326
      }
2278 2327
      if (_node_potential) {
2279 2328
        delete _node_potential;
2280 2329
      }
2281 2330
      if (_blossom_set) {
2282 2331
        delete _blossom_index;
2283 2332
        delete _blossom_set;
2284 2333
        delete _blossom_data;
2285 2334
      }
2286 2335

	
2287 2336
      if (_node_index) {
2288 2337
        delete _node_index;
2289 2338
        delete _node_heap_index;
2290 2339
        delete _node_data;
2291 2340
      }
2292 2341

	
2293 2342
      if (_tree_set) {
2294 2343
        delete _tree_set_index;
2295 2344
        delete _tree_set;
2296 2345
      }
2297 2346
      if (_delta2) {
2298 2347
        delete _delta2_index;
2299 2348
        delete _delta2;
2300 2349
      }
2301 2350
      if (_delta3) {
2302 2351
        delete _delta3_index;
2303 2352
        delete _delta3;
2304 2353
      }
2305 2354
      if (_delta4) {
2306 2355
        delete _delta4_index;
2307 2356
        delete _delta4;
2308 2357
      }
2309 2358
    }
2310 2359

	
2311 2360
    void matchedToEven(int blossom, int tree) {
2312 2361
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2313 2362
        _delta2->erase(blossom);
2314 2363
      }
2315 2364

	
2316 2365
      if (!_blossom_set->trivial(blossom)) {
2317 2366
        (*_blossom_data)[blossom].pot -=
2318 2367
          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2319 2368
      }
2320 2369

	
2321 2370
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2322 2371
           n != INVALID; ++n) {
2323 2372

	
2324 2373
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
2325 2374
        int ni = (*_node_index)[n];
2326 2375

	
2327 2376
        (*_node_data)[ni].heap.clear();
2328 2377
        (*_node_data)[ni].heap_index.clear();
2329 2378

	
2330 2379
        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
2331 2380

	
2332 2381
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
2333 2382
          Node v = _graph.source(e);
2334 2383
          int vb = _blossom_set->find(v);
2335 2384
          int vi = (*_node_index)[v];
2336 2385

	
2337 2386
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2338 2387
            dualScale * _weight[e];
2339 2388

	
2340 2389
          if ((*_blossom_data)[vb].status == EVEN) {
2341 2390
            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2342 2391
              _delta3->push(e, rw / 2);
2343 2392
            }
2344 2393
          } else {
2345 2394
            typename std::map<int, Arc>::iterator it =
2346 2395
              (*_node_data)[vi].heap_index.find(tree);
2347 2396

	
2348 2397
            if (it != (*_node_data)[vi].heap_index.end()) {
2349 2398
              if ((*_node_data)[vi].heap[it->second] > rw) {
2350 2399
                (*_node_data)[vi].heap.replace(it->second, e);
2351 2400
                (*_node_data)[vi].heap.decrease(e, rw);
2352 2401
                it->second = e;
2353 2402
              }
2354 2403
            } else {
2355 2404
              (*_node_data)[vi].heap.push(e, rw);
2356 2405
              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2357 2406
            }
2358 2407

	
2359 2408
            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2360 2409
              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2361 2410

	
2362 2411
              if ((*_blossom_data)[vb].status == MATCHED) {
2363 2412
                if (_delta2->state(vb) != _delta2->IN_HEAP) {
2364 2413
                  _delta2->push(vb, _blossom_set->classPrio(vb) -
2365 2414
                               (*_blossom_data)[vb].offset);
2366 2415
                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2367 2416
                           (*_blossom_data)[vb].offset){
2368 2417
                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2369 2418
                                   (*_blossom_data)[vb].offset);
2370 2419
                }
2371 2420
              }
2372 2421
            }
2373 2422
          }
2374 2423
        }
2375 2424
      }
2376 2425
      (*_blossom_data)[blossom].offset = 0;
2377 2426
    }
2378 2427

	
2379 2428
    void matchedToOdd(int blossom) {
2380 2429
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2381 2430
        _delta2->erase(blossom);
2382 2431
      }
2383 2432
      (*_blossom_data)[blossom].offset += _delta_sum;
2384 2433
      if (!_blossom_set->trivial(blossom)) {
2385 2434
        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
2386 2435
                     (*_blossom_data)[blossom].offset);
2387 2436
      }
2388 2437
    }
2389 2438

	
2390 2439
    void evenToMatched(int blossom, int tree) {
2391 2440
      if (!_blossom_set->trivial(blossom)) {
2392 2441
        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
2393 2442
      }
2394 2443

	
2395 2444
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2396 2445
           n != INVALID; ++n) {
2397 2446
        int ni = (*_node_index)[n];
2398 2447
        (*_node_data)[ni].pot -= _delta_sum;
2399 2448

	
2400 2449
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
2401 2450
          Node v = _graph.source(e);
2402 2451
          int vb = _blossom_set->find(v);
2403 2452
          int vi = (*_node_index)[v];
2404 2453

	
2405 2454
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2406 2455
            dualScale * _weight[e];
2407 2456

	
2408 2457
          if (vb == blossom) {
2409 2458
            if (_delta3->state(e) == _delta3->IN_HEAP) {
2410 2459
              _delta3->erase(e);
2411 2460
            }
2412 2461
          } else if ((*_blossom_data)[vb].status == EVEN) {
2413 2462

	
2414 2463
            if (_delta3->state(e) == _delta3->IN_HEAP) {
2415 2464
              _delta3->erase(e);
2416 2465
            }
2417 2466

	
2418 2467
            int vt = _tree_set->find(vb);
2419 2468

	
2420 2469
            if (vt != tree) {
2421 2470

	
2422 2471
              Arc r = _graph.oppositeArc(e);
2423 2472

	
2424 2473
              typename std::map<int, Arc>::iterator it =
2425 2474
                (*_node_data)[ni].heap_index.find(vt);
2426 2475

	
2427 2476
              if (it != (*_node_data)[ni].heap_index.end()) {
2428 2477
                if ((*_node_data)[ni].heap[it->second] > rw) {
2429 2478
                  (*_node_data)[ni].heap.replace(it->second, r);
2430 2479
                  (*_node_data)[ni].heap.decrease(r, rw);
2431 2480
                  it->second = r;
2432 2481
                }
2433 2482
              } else {
2434 2483
                (*_node_data)[ni].heap.push(r, rw);
2435 2484
                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2436 2485
              }
2437 2486

	
2438 2487
              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2439 2488
                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2440 2489

	
2441 2490
                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2442 2491
                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2443 2492
                               (*_blossom_data)[blossom].offset);
2444 2493
                } else if ((*_delta2)[blossom] >
2445 2494
                           _blossom_set->classPrio(blossom) -
2446 2495
                           (*_blossom_data)[blossom].offset){
2447 2496
                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2448 2497
                                   (*_blossom_data)[blossom].offset);
2449 2498
                }
2450 2499
              }
2451 2500
            }
2452 2501
          } else {
2453 2502

	
2454 2503
            typename std::map<int, Arc>::iterator it =
2455 2504
              (*_node_data)[vi].heap_index.find(tree);
2456 2505

	
2457 2506
            if (it != (*_node_data)[vi].heap_index.end()) {
2458 2507
              (*_node_data)[vi].heap.erase(it->second);
2459 2508
              (*_node_data)[vi].heap_index.erase(it);
2460 2509
              if ((*_node_data)[vi].heap.empty()) {
2461 2510
                _blossom_set->increase(v, std::numeric_limits<Value>::max());
2462 2511
              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
2463 2512
                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
2464 2513
              }
2465 2514

	
2466 2515
              if ((*_blossom_data)[vb].status == MATCHED) {
2467 2516
                if (_blossom_set->classPrio(vb) ==
2468 2517
                    std::numeric_limits<Value>::max()) {
2469 2518
                  _delta2->erase(vb);
2470 2519
                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
2471 2520
                           (*_blossom_data)[vb].offset) {
2472 2521
                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
2473 2522
                                   (*_blossom_data)[vb].offset);
2474 2523
                }
2475 2524
              }
2476 2525
            }
2477 2526
          }
2478 2527
        }
2479 2528
      }
2480 2529
    }
2481 2530

	
2482 2531
    void oddToMatched(int blossom) {
2483 2532
      (*_blossom_data)[blossom].offset -= _delta_sum;
2484 2533

	
2485 2534
      if (_blossom_set->classPrio(blossom) !=
2486 2535
          std::numeric_limits<Value>::max()) {
2487 2536
        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2488 2537
                       (*_blossom_data)[blossom].offset);
2489 2538
      }
2490 2539

	
2491 2540
      if (!_blossom_set->trivial(blossom)) {
2492 2541
        _delta4->erase(blossom);
2493 2542
      }
2494 2543
    }
2495 2544

	
2496 2545
    void oddToEven(int blossom, int tree) {
2497 2546
      if (!_blossom_set->trivial(blossom)) {
2498 2547
        _delta4->erase(blossom);
2499 2548
        (*_blossom_data)[blossom].pot -=
2500 2549
          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
2501 2550
      }
2502 2551

	
2503 2552
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2504 2553
           n != INVALID; ++n) {
2505 2554
        int ni = (*_node_index)[n];
2506 2555

	
2507 2556
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
2508 2557

	
2509 2558
        (*_node_data)[ni].heap.clear();
2510 2559
        (*_node_data)[ni].heap_index.clear();
2511 2560
        (*_node_data)[ni].pot +=
2512 2561
          2 * _delta_sum - (*_blossom_data)[blossom].offset;
2513 2562

	
2514 2563
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
2515 2564
          Node v = _graph.source(e);
2516 2565
          int vb = _blossom_set->find(v);
2517 2566
          int vi = (*_node_index)[v];
2518 2567

	
2519 2568
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2520 2569
            dualScale * _weight[e];
2521 2570

	
2522 2571
          if ((*_blossom_data)[vb].status == EVEN) {
2523 2572
            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2524 2573
              _delta3->push(e, rw / 2);
2525 2574
            }
2526 2575
          } else {
2527 2576

	
2528 2577
            typename std::map<int, Arc>::iterator it =
2529 2578
              (*_node_data)[vi].heap_index.find(tree);
2530 2579

	
2531 2580
            if (it != (*_node_data)[vi].heap_index.end()) {
2532 2581
              if ((*_node_data)[vi].heap[it->second] > rw) {
2533 2582
                (*_node_data)[vi].heap.replace(it->second, e);
2534 2583
                (*_node_data)[vi].heap.decrease(e, rw);
2535 2584
                it->second = e;
2536 2585
              }
2537 2586
            } else {
2538 2587
              (*_node_data)[vi].heap.push(e, rw);
2539 2588
              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2540 2589
            }
2541 2590

	
2542 2591
            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2543 2592
              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2544 2593

	
2545 2594
              if ((*_blossom_data)[vb].status == MATCHED) {
2546 2595
                if (_delta2->state(vb) != _delta2->IN_HEAP) {
2547 2596
                  _delta2->push(vb, _blossom_set->classPrio(vb) -
2548 2597
                               (*_blossom_data)[vb].offset);
2549 2598
                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2550 2599
                           (*_blossom_data)[vb].offset) {
2551 2600
                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2552 2601
                                   (*_blossom_data)[vb].offset);
2553 2602
                }
2554 2603
              }
2555 2604
            }
2556 2605
          }
2557 2606
        }
2558 2607
      }
2559 2608
      (*_blossom_data)[blossom].offset = 0;
2560 2609
    }
2561 2610

	
2562 2611
    void alternatePath(int even, int tree) {
2563 2612
      int odd;
2564 2613

	
2565 2614
      evenToMatched(even, tree);
2566 2615
      (*_blossom_data)[even].status = MATCHED;
2567 2616

	
2568 2617
      while ((*_blossom_data)[even].pred != INVALID) {
2569 2618
        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
2570 2619
        (*_blossom_data)[odd].status = MATCHED;
2571 2620
        oddToMatched(odd);
2572 2621
        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2573 2622

	
2574 2623
        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
2575 2624
        (*_blossom_data)[even].status = MATCHED;
2576 2625
        evenToMatched(even, tree);
2577 2626
        (*_blossom_data)[even].next =
2578 2627
          _graph.oppositeArc((*_blossom_data)[odd].pred);
2579 2628
      }
2580 2629

	
2581 2630
    }
2582 2631

	
2583 2632
    void destroyTree(int tree) {
2584 2633
      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2585 2634
        if ((*_blossom_data)[b].status == EVEN) {
2586 2635
          (*_blossom_data)[b].status = MATCHED;
2587 2636
          evenToMatched(b, tree);
2588 2637
        } else if ((*_blossom_data)[b].status == ODD) {
2589 2638
          (*_blossom_data)[b].status = MATCHED;
2590 2639
          oddToMatched(b);
2591 2640
        }
2592 2641
      }
2593 2642
      _tree_set->eraseClass(tree);
2594 2643
    }
2595 2644

	
2596 2645
    void augmentOnEdge(const Edge& edge) {
2597 2646

	
2598 2647
      int left = _blossom_set->find(_graph.u(edge));
2599 2648
      int right = _blossom_set->find(_graph.v(edge));
2600 2649

	
2601 2650
      int left_tree = _tree_set->find(left);
2602 2651
      alternatePath(left, left_tree);
2603 2652
      destroyTree(left_tree);
2604 2653

	
2605 2654
      int right_tree = _tree_set->find(right);
2606 2655
      alternatePath(right, right_tree);
2607 2656
      destroyTree(right_tree);
2608 2657

	
2609 2658
      (*_blossom_data)[left].next = _graph.direct(edge, true);
2610 2659
      (*_blossom_data)[right].next = _graph.direct(edge, false);
2611 2660
    }
2612 2661

	
2613 2662
    void extendOnArc(const Arc& arc) {
2614 2663
      int base = _blossom_set->find(_graph.target(arc));
2615 2664
      int tree = _tree_set->find(base);
2616 2665

	
2617 2666
      int odd = _blossom_set->find(_graph.source(arc));
2618 2667
      _tree_set->insert(odd, tree);
2619 2668
      (*_blossom_data)[odd].status = ODD;
2620 2669
      matchedToOdd(odd);
2621 2670
      (*_blossom_data)[odd].pred = arc;
2622 2671

	
2623 2672
      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
2624 2673
      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2625 2674
      _tree_set->insert(even, tree);
2626 2675
      (*_blossom_data)[even].status = EVEN;
2627 2676
      matchedToEven(even, tree);
2628 2677
    }
2629 2678

	
2630 2679
    void shrinkOnEdge(const Edge& edge, int tree) {
2631 2680
      int nca = -1;
2632 2681
      std::vector<int> left_path, right_path;
2633 2682

	
2634 2683
      {
2635 2684
        std::set<int> left_set, right_set;
2636 2685
        int left = _blossom_set->find(_graph.u(edge));
2637 2686
        left_path.push_back(left);
2638 2687
        left_set.insert(left);
2639 2688

	
2640 2689
        int right = _blossom_set->find(_graph.v(edge));
2641 2690
        right_path.push_back(right);
2642 2691
        right_set.insert(right);
2643 2692

	
2644 2693
        while (true) {
2645 2694

	
2646 2695
          if ((*_blossom_data)[left].pred == INVALID) break;
2647 2696

	
2648 2697
          left =
2649 2698
            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2650 2699
          left_path.push_back(left);
2651 2700
          left =
2652 2701
            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2653 2702
          left_path.push_back(left);
2654 2703

	
2655 2704
          left_set.insert(left);
2656 2705

	
2657 2706
          if (right_set.find(left) != right_set.end()) {
2658 2707
            nca = left;
2659 2708
            break;
2660 2709
          }
2661 2710

	
2662 2711
          if ((*_blossom_data)[right].pred == INVALID) break;
2663 2712

	
2664 2713
          right =
2665 2714
            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2666 2715
          right_path.push_back(right);
2667 2716
          right =
2668 2717
            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2669 2718
          right_path.push_back(right);
2670 2719

	
2671 2720
          right_set.insert(right);
2672 2721

	
2673 2722
          if (left_set.find(right) != left_set.end()) {
2674 2723
            nca = right;
2675 2724
            break;
2676 2725
          }
2677 2726

	
2678 2727
        }
2679 2728

	
2680 2729
        if (nca == -1) {
2681 2730
          if ((*_blossom_data)[left].pred == INVALID) {
2682 2731
            nca = right;
2683 2732
            while (left_set.find(nca) == left_set.end()) {
2684 2733
              nca =
2685 2734
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2686 2735
              right_path.push_back(nca);
2687 2736
              nca =
2688 2737
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2689 2738
              right_path.push_back(nca);
2690 2739
            }
2691 2740
          } else {
2692 2741
            nca = left;
2693 2742
            while (right_set.find(nca) == right_set.end()) {
2694 2743
              nca =
2695 2744
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2696 2745
              left_path.push_back(nca);
2697 2746
              nca =
2698 2747
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2699 2748
              left_path.push_back(nca);
2700 2749
            }
2701 2750
          }
2702 2751
        }
2703 2752
      }
2704 2753

	
2705 2754
      std::vector<int> subblossoms;
2706 2755
      Arc prev;
2707 2756

	
2708 2757
      prev = _graph.direct(edge, true);
2709 2758
      for (int i = 0; left_path[i] != nca; i += 2) {
2710 2759
        subblossoms.push_back(left_path[i]);
2711 2760
        (*_blossom_data)[left_path[i]].next = prev;
2712 2761
        _tree_set->erase(left_path[i]);
2713 2762

	
2714 2763
        subblossoms.push_back(left_path[i + 1]);
2715 2764
        (*_blossom_data)[left_path[i + 1]].status = EVEN;
2716 2765
        oddToEven(left_path[i + 1], tree);
2717 2766
        _tree_set->erase(left_path[i + 1]);
2718 2767
        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
2719 2768
      }
2720 2769

	
2721 2770
      int k = 0;
2722 2771
      while (right_path[k] != nca) ++k;
2723 2772

	
2724 2773
      subblossoms.push_back(nca);
2725 2774
      (*_blossom_data)[nca].next = prev;
2726 2775

	
2727 2776
      for (int i = k - 2; i >= 0; i -= 2) {
2728 2777
        subblossoms.push_back(right_path[i + 1]);
2729 2778
        (*_blossom_data)[right_path[i + 1]].status = EVEN;
2730 2779
        oddToEven(right_path[i + 1], tree);
2731 2780
        _tree_set->erase(right_path[i + 1]);
2732 2781

	
2733 2782
        (*_blossom_data)[right_path[i + 1]].next =
2734 2783
          (*_blossom_data)[right_path[i + 1]].pred;
2735 2784

	
2736 2785
        subblossoms.push_back(right_path[i]);
2737 2786
        _tree_set->erase(right_path[i]);
2738 2787
      }
2739 2788

	
2740 2789
      int surface =
2741 2790
        _blossom_set->join(subblossoms.begin(), subblossoms.end());
2742 2791

	
2743 2792
      for (int i = 0; i < int(subblossoms.size()); ++i) {
2744 2793
        if (!_blossom_set->trivial(subblossoms[i])) {
2745 2794
          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
2746 2795
        }
2747 2796
        (*_blossom_data)[subblossoms[i]].status = MATCHED;
2748 2797
      }
2749 2798

	
2750 2799
      (*_blossom_data)[surface].pot = -2 * _delta_sum;
2751 2800
      (*_blossom_data)[surface].offset = 0;
2752 2801
      (*_blossom_data)[surface].status = EVEN;
2753 2802
      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
2754 2803
      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
2755 2804

	
2756 2805
      _tree_set->insert(surface, tree);
2757 2806
      _tree_set->erase(nca);
2758 2807
    }
2759 2808

	
2760 2809
    void splitBlossom(int blossom) {
2761 2810
      Arc next = (*_blossom_data)[blossom].next;
2762 2811
      Arc pred = (*_blossom_data)[blossom].pred;
2763 2812

	
2764 2813
      int tree = _tree_set->find(blossom);
2765 2814

	
2766 2815
      (*_blossom_data)[blossom].status = MATCHED;
2767 2816
      oddToMatched(blossom);
2768 2817
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2769 2818
        _delta2->erase(blossom);
2770 2819
      }
2771 2820

	
2772 2821
      std::vector<int> subblossoms;
2773 2822
      _blossom_set->split(blossom, std::back_inserter(subblossoms));
2774 2823

	
2775 2824
      Value offset = (*_blossom_data)[blossom].offset;
2776 2825
      int b = _blossom_set->find(_graph.source(pred));
2777 2826
      int d = _blossom_set->find(_graph.source(next));
2778 2827

	
2779 2828
      int ib = -1, id = -1;
2780 2829
      for (int i = 0; i < int(subblossoms.size()); ++i) {
2781 2830
        if (subblossoms[i] == b) ib = i;
2782 2831
        if (subblossoms[i] == d) id = i;
2783 2832

	
2784 2833
        (*_blossom_data)[subblossoms[i]].offset = offset;
2785 2834
        if (!_blossom_set->trivial(subblossoms[i])) {
2786 2835
          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2787 2836
        }
2788 2837
        if (_blossom_set->classPrio(subblossoms[i]) !=
2789 2838
            std::numeric_limits<Value>::max()) {
2790 2839
          _delta2->push(subblossoms[i],
2791 2840
                        _blossom_set->classPrio(subblossoms[i]) -
2792 2841
                        (*_blossom_data)[subblossoms[i]].offset);
2793 2842
        }
2794 2843
      }
2795 2844

	
2796 2845
      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2797 2846
        for (int i = (id + 1) % subblossoms.size();
2798 2847
             i != ib; i = (i + 2) % subblossoms.size()) {
2799 2848
          int sb = subblossoms[i];
2800 2849
          int tb = subblossoms[(i + 1) % subblossoms.size()];
2801 2850
          (*_blossom_data)[sb].next =
2802 2851
            _graph.oppositeArc((*_blossom_data)[tb].next);
2803 2852
        }
2804 2853

	
2805 2854
        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2806 2855
          int sb = subblossoms[i];
2807 2856
          int tb = subblossoms[(i + 1) % subblossoms.size()];
2808 2857
          int ub = subblossoms[(i + 2) % subblossoms.size()];
2809 2858

	
2810 2859
          (*_blossom_data)[sb].status = ODD;
2811 2860
          matchedToOdd(sb);
2812 2861
          _tree_set->insert(sb, tree);
2813 2862
          (*_blossom_data)[sb].pred = pred;
2814 2863
          (*_blossom_data)[sb].next =
2815 2864
                           _graph.oppositeArc((*_blossom_data)[tb].next);
2816 2865

	
2817 2866
          pred = (*_blossom_data)[ub].next;
2818 2867

	
2819 2868
          (*_blossom_data)[tb].status = EVEN;
2820 2869
          matchedToEven(tb, tree);
2821 2870
          _tree_set->insert(tb, tree);
2822 2871
          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2823 2872
        }
2824 2873

	
2825 2874
        (*_blossom_data)[subblossoms[id]].status = ODD;
2826 2875
        matchedToOdd(subblossoms[id]);
2827 2876
        _tree_set->insert(subblossoms[id], tree);
2828 2877
        (*_blossom_data)[subblossoms[id]].next = next;
2829 2878
        (*_blossom_data)[subblossoms[id]].pred = pred;
2830 2879

	
2831 2880
      } else {
2832 2881

	
2833 2882
        for (int i = (ib + 1) % subblossoms.size();
2834 2883
             i != id; i = (i + 2) % subblossoms.size()) {
2835 2884
          int sb = subblossoms[i];
2836 2885
          int tb = subblossoms[(i + 1) % subblossoms.size()];
2837 2886
          (*_blossom_data)[sb].next =
2838 2887
            _graph.oppositeArc((*_blossom_data)[tb].next);
2839 2888
        }
2840 2889

	
2841 2890
        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2842 2891
          int sb = subblossoms[i];
2843 2892
          int tb = subblossoms[(i + 1) % subblossoms.size()];
2844 2893
          int ub = subblossoms[(i + 2) % subblossoms.size()];
2845 2894

	
2846 2895
          (*_blossom_data)[sb].status = ODD;
2847 2896
          matchedToOdd(sb);
2848 2897
          _tree_set->insert(sb, tree);
2849 2898
          (*_blossom_data)[sb].next = next;
2850 2899
          (*_blossom_data)[sb].pred =
2851 2900
            _graph.oppositeArc((*_blossom_data)[tb].next);
2852 2901

	
2853 2902
          (*_blossom_data)[tb].status = EVEN;
2854 2903
          matchedToEven(tb, tree);
2855 2904
          _tree_set->insert(tb, tree);
2856 2905
          (*_blossom_data)[tb].pred =
2857 2906
            (*_blossom_data)[tb].next =
2858 2907
            _graph.oppositeArc((*_blossom_data)[ub].next);
2859 2908
          next = (*_blossom_data)[ub].next;
2860 2909
        }
2861 2910

	
2862 2911
        (*_blossom_data)[subblossoms[ib]].status = ODD;
2863 2912
        matchedToOdd(subblossoms[ib]);
2864 2913
        _tree_set->insert(subblossoms[ib], tree);
2865 2914
        (*_blossom_data)[subblossoms[ib]].next = next;
2866 2915
        (*_blossom_data)[subblossoms[ib]].pred = pred;
2867 2916
      }
2868 2917
      _tree_set->erase(blossom);
2869 2918
    }
2870 2919

	
2871 2920
    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2872 2921
      if (_blossom_set->trivial(blossom)) {
2873 2922
        int bi = (*_node_index)[base];
2874 2923
        Value pot = (*_node_data)[bi].pot;
2875 2924

	
2876 2925
        (*_matching)[base] = matching;
2877 2926
        _blossom_node_list.push_back(base);
2878 2927
        (*_node_potential)[base] = pot;
2879 2928
      } else {
2880 2929

	
2881 2930
        Value pot = (*_blossom_data)[blossom].pot;
2882 2931
        int bn = _blossom_node_list.size();
2883 2932

	
2884 2933
        std::vector<int> subblossoms;
2885 2934
        _blossom_set->split(blossom, std::back_inserter(subblossoms));
2886 2935
        int b = _blossom_set->find(base);
2887 2936
        int ib = -1;
2888 2937
        for (int i = 0; i < int(subblossoms.size()); ++i) {
2889 2938
          if (subblossoms[i] == b) { ib = i; break; }
2890 2939
        }
2891 2940

	
2892 2941
        for (int i = 1; i < int(subblossoms.size()); i += 2) {
2893 2942
          int sb = subblossoms[(ib + i) % subblossoms.size()];
2894 2943
          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2895 2944

	
2896 2945
          Arc m = (*_blossom_data)[tb].next;
2897 2946
          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2898 2947
          extractBlossom(tb, _graph.source(m), m);
2899 2948
        }
2900 2949
        extractBlossom(subblossoms[ib], base, matching);
2901 2950

	
2902 2951
        int en = _blossom_node_list.size();
2903 2952

	
2904 2953
        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2905 2954
      }
2906 2955
    }
2907 2956

	
2908 2957
    void extractMatching() {
2909 2958
      std::vector<int> blossoms;
2910 2959
      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2911 2960
        blossoms.push_back(c);
2912 2961
      }
2913 2962

	
2914 2963
      for (int i = 0; i < int(blossoms.size()); ++i) {
2915 2964

	
2916 2965
        Value offset = (*_blossom_data)[blossoms[i]].offset;
2917 2966
        (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2918 2967
        for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2919 2968
             n != INVALID; ++n) {
2920 2969
          (*_node_data)[(*_node_index)[n]].pot -= offset;
2921 2970
        }
2922 2971

	
2923 2972
        Arc matching = (*_blossom_data)[blossoms[i]].next;
2924 2973
        Node base = _graph.source(matching);
2925 2974
        extractBlossom(blossoms[i], base, matching);
2926 2975
      }
2927 2976
    }
2928 2977

	
2929 2978
  public:
2930 2979

	
2931 2980
    /// \brief Constructor
2932 2981
    ///
2933 2982
    /// Constructor.
2934 2983
    MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2935 2984
      : _graph(graph), _weight(weight), _matching(0),
2936 2985
        _node_potential(0), _blossom_potential(), _blossom_node_list(),
2937 2986
        _node_num(0), _blossom_num(0),
2938 2987

	
2939 2988
        _blossom_index(0), _blossom_set(0), _blossom_data(0),
2940 2989
        _node_index(0), _node_heap_index(0), _node_data(0),
2941 2990
        _tree_set_index(0), _tree_set(0),
2942 2991

	
2943 2992
        _delta2_index(0), _delta2(0),
2944 2993
        _delta3_index(0), _delta3(0),
2945 2994
        _delta4_index(0), _delta4(0),
2946 2995

	
2947 2996
        _delta_sum(), _unmatched(0),
2948 2997

	
2949 2998
        _fractional(0)
2950 2999
    {}
2951 3000

	
2952 3001
    ~MaxWeightedPerfectMatching() {
2953 3002
      destroyStructures();
2954 3003
      if (_fractional) {
2955 3004
        delete _fractional;
2956 3005
      }
2957 3006
    }
2958 3007

	
2959 3008
    /// \name Execution Control
2960 3009
    /// The simplest way to execute the algorithm is to use the
2961 3010
    /// \ref run() member function.
2962 3011

	
2963 3012
    ///@{
2964 3013

	
2965 3014
    /// \brief Initialize the algorithm
2966 3015
    ///
2967 3016
    /// This function initializes the algorithm.
2968 3017
    void init() {
2969 3018
      createStructures();
2970 3019

	
3020
      _blossom_node_list.clear();
3021
      _blossom_potential.clear();
3022

	
2971 3023
      for (ArcIt e(_graph); e != INVALID; ++e) {
2972 3024
        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
2973 3025
      }
2974 3026
      for (EdgeIt e(_graph); e != INVALID; ++e) {
2975 3027
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
2976 3028
      }
2977 3029
      for (int i = 0; i < _blossom_num; ++i) {
2978 3030
        (*_delta2_index)[i] = _delta2->PRE_HEAP;
2979 3031
        (*_delta4_index)[i] = _delta4->PRE_HEAP;
2980 3032
      }
2981 3033

	
2982 3034
      _unmatched = _node_num;
2983 3035

	
3036
      _delta2->clear();
3037
      _delta3->clear();
3038
      _delta4->clear();
3039
      _blossom_set->clear();
3040
      _tree_set->clear();
3041

	
2984 3042
      int index = 0;
2985 3043
      for (NodeIt n(_graph); n != INVALID; ++n) {
2986 3044
        Value max = - std::numeric_limits<Value>::max();
2987 3045
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2988 3046
          if (_graph.target(e) == n) continue;
2989 3047
          if ((dualScale * _weight[e]) / 2 > max) {
2990 3048
            max = (dualScale * _weight[e]) / 2;
2991 3049
          }
2992 3050
        }
2993 3051
        (*_node_index)[n] = index;
3052
        (*_node_data)[index].heap_index.clear();
3053
        (*_node_data)[index].heap.clear();
2994 3054
        (*_node_data)[index].pot = max;
2995 3055
        int blossom =
2996 3056
          _blossom_set->insert(n, std::numeric_limits<Value>::max());
2997 3057

	
2998 3058
        _tree_set->insert(blossom);
2999 3059

	
3000 3060
        (*_blossom_data)[blossom].status = EVEN;
3001 3061
        (*_blossom_data)[blossom].pred = INVALID;
3002 3062
        (*_blossom_data)[blossom].next = INVALID;
3003 3063
        (*_blossom_data)[blossom].pot = 0;
3004 3064
        (*_blossom_data)[blossom].offset = 0;
3005 3065
        ++index;
3006 3066
      }
3007 3067
      for (EdgeIt e(_graph); e != INVALID; ++e) {
3008 3068
        int si = (*_node_index)[_graph.u(e)];
3009 3069
        int ti = (*_node_index)[_graph.v(e)];
3010 3070
        if (_graph.u(e) != _graph.v(e)) {
3011 3071
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
3012 3072
                            dualScale * _weight[e]) / 2);
3013 3073
        }
3014 3074
      }
3015 3075
    }
3016 3076

	
3017 3077
    /// \brief Initialize the algorithm with fractional matching
3018 3078
    ///
3019 3079
    /// This function initializes the algorithm with a fractional
3020 3080
    /// matching. This initialization is also called jumpstart heuristic.
3021 3081
    void fractionalInit() {
3022 3082
      createStructures();
3023 3083
      
3024 3084
      if (_fractional == 0) {
3025 3085
        _fractional = new FractionalMatching(_graph, _weight, false);
3026 3086
      }
3027 3087
      if (!_fractional->run()) {
3028 3088
        _unmatched = -1;
3029 3089
        return;
3030 3090
      }
3031 3091

	
3032 3092
      for (ArcIt e(_graph); e != INVALID; ++e) {
3033 3093
        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
3034 3094
      }
3035 3095
      for (EdgeIt e(_graph); e != INVALID; ++e) {
3036 3096
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
3037 3097
      }
3038 3098
      for (int i = 0; i < _blossom_num; ++i) {
3039 3099
        (*_delta2_index)[i] = _delta2->PRE_HEAP;
3040 3100
        (*_delta4_index)[i] = _delta4->PRE_HEAP;
3041 3101
      }
3042 3102

	
3043 3103
      _unmatched = 0;
3044 3104

	
3045 3105
      int index = 0;
3046 3106
      for (NodeIt n(_graph); n != INVALID; ++n) {
3047 3107
        Value pot = _fractional->nodeValue(n);
3048 3108
        (*_node_index)[n] = index;
3049 3109
        (*_node_data)[index].pot = pot;
3050 3110
        int blossom =
3051 3111
          _blossom_set->insert(n, std::numeric_limits<Value>::max());
3052 3112

	
3053 3113
        (*_blossom_data)[blossom].status = MATCHED;
3054 3114
        (*_blossom_data)[blossom].pred = INVALID;
3055 3115
        (*_blossom_data)[blossom].next = _fractional->matching(n);
3056 3116
        (*_blossom_data)[blossom].pot = 0;
3057 3117
        (*_blossom_data)[blossom].offset = 0;
3058 3118
        ++index;
3059 3119
      }
3060 3120

	
3061 3121
      typename Graph::template NodeMap<bool> processed(_graph, false);
3062 3122
      for (NodeIt n(_graph); n != INVALID; ++n) {
3063 3123
        if (processed[n]) continue;
3064 3124
        processed[n] = true;
3065 3125
        if (_fractional->matching(n) == INVALID) continue;
3066 3126
        int num = 1;
3067 3127
        Node v = _graph.target(_fractional->matching(n));
3068 3128
        while (n != v) {
3069 3129
          processed[v] = true;
3070 3130
          v = _graph.target(_fractional->matching(v));
3071 3131
          ++num;
3072 3132
        }
3073 3133

	
3074 3134
        if (num % 2 == 1) {
3075 3135
          std::vector<int> subblossoms(num);
3076 3136

	
3077 3137
          subblossoms[--num] = _blossom_set->find(n);
3078 3138
          v = _graph.target(_fractional->matching(n));
3079 3139
          while (n != v) {
3080 3140
            subblossoms[--num] = _blossom_set->find(v);
3081 3141
            v = _graph.target(_fractional->matching(v));            
3082 3142
          }
3083 3143
          
3084 3144
          int surface = 
3085 3145
            _blossom_set->join(subblossoms.begin(), subblossoms.end());
3086 3146
          (*_blossom_data)[surface].status = EVEN;
3087 3147
          (*_blossom_data)[surface].pred = INVALID;
3088 3148
          (*_blossom_data)[surface].next = INVALID;
3089 3149
          (*_blossom_data)[surface].pot = 0;
3090 3150
          (*_blossom_data)[surface].offset = 0;
3091 3151
          
3092 3152
          _tree_set->insert(surface);
3093 3153
          ++_unmatched;
3094 3154
        }
3095 3155
      }
3096 3156

	
3097 3157
      for (EdgeIt e(_graph); e != INVALID; ++e) {
3098 3158
        int si = (*_node_index)[_graph.u(e)];
3099 3159
        int sb = _blossom_set->find(_graph.u(e));
3100 3160
        int ti = (*_node_index)[_graph.v(e)];
3101 3161
        int tb = _blossom_set->find(_graph.v(e));
3102 3162
        if ((*_blossom_data)[sb].status == EVEN &&
3103 3163
            (*_blossom_data)[tb].status == EVEN && sb != tb) {
3104 3164
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
3105 3165
                            dualScale * _weight[e]) / 2);
3106 3166
        }
3107 3167
      }
3108 3168

	
3109 3169
      for (NodeIt n(_graph); n != INVALID; ++n) {
3110 3170
        int nb = _blossom_set->find(n);
3111 3171
        if ((*_blossom_data)[nb].status != MATCHED) continue;
3112 3172
        int ni = (*_node_index)[n];
3113 3173

	
3114 3174
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
3115 3175
          Node v = _graph.target(e);
3116 3176
          int vb = _blossom_set->find(v);
3117 3177
          int vi = (*_node_index)[v];
3118 3178

	
3119 3179
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
3120 3180
            dualScale * _weight[e];
3121 3181

	
3122 3182
          if ((*_blossom_data)[vb].status == EVEN) {
3123 3183

	
3124 3184
            int vt = _tree_set->find(vb);
3125 3185

	
3126 3186
            typename std::map<int, Arc>::iterator it =
3127 3187
              (*_node_data)[ni].heap_index.find(vt);
3128 3188

	
3129 3189
            if (it != (*_node_data)[ni].heap_index.end()) {
3130 3190
              if ((*_node_data)[ni].heap[it->second] > rw) {
3131 3191
                (*_node_data)[ni].heap.replace(it->second, e);
3132 3192
                (*_node_data)[ni].heap.decrease(e, rw);
3133 3193
                it->second = e;
3134 3194
              }
3135 3195
            } else {
3136 3196
              (*_node_data)[ni].heap.push(e, rw);
3137 3197
              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
3138 3198
            }
3139 3199
          }
3140 3200
        }
3141 3201
            
3142 3202
        if (!(*_node_data)[ni].heap.empty()) {
3143 3203
          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
3144 3204
          _delta2->push(nb, _blossom_set->classPrio(nb));
3145 3205
        }
3146 3206
      }
3147 3207
    }
3148 3208

	
3149 3209
    /// \brief Start the algorithm
3150 3210
    ///
3151 3211
    /// This function starts the algorithm.
3152 3212
    ///
3153 3213
    /// \pre \ref init() or \ref fractionalInit() must be called before
3154 3214
    /// using this function.
3155 3215
    bool start() {
3156 3216
      enum OpType {
3157 3217
        D2, D3, D4
3158 3218
      };
3159 3219

	
3160 3220
      if (_unmatched == -1) return false;
3161 3221

	
3162 3222
      while (_unmatched > 0) {
3163 3223
        Value d2 = !_delta2->empty() ?
3164 3224
          _delta2->prio() : std::numeric_limits<Value>::max();
3165 3225

	
3166 3226
        Value d3 = !_delta3->empty() ?
3167 3227
          _delta3->prio() : std::numeric_limits<Value>::max();
3168 3228

	
3169 3229
        Value d4 = !_delta4->empty() ?
3170 3230
          _delta4->prio() : std::numeric_limits<Value>::max();
3171 3231

	
3172 3232
        _delta_sum = d3; OpType ot = D3;
3173 3233
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
3174 3234
        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
3175 3235

	
3176 3236
        if (_delta_sum == std::numeric_limits<Value>::max()) {
3177 3237
          return false;
3178 3238
        }
3179 3239

	
3180 3240
        switch (ot) {
3181 3241
        case D2:
3182 3242
          {
3183 3243
            int blossom = _delta2->top();
3184 3244
            Node n = _blossom_set->classTop(blossom);
3185 3245
            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
3186 3246
            extendOnArc(e);
3187 3247
          }
3188 3248
          break;
3189 3249
        case D3:
3190 3250
          {
3191 3251
            Edge e = _delta3->top();
3192 3252

	
3193 3253
            int left_blossom = _blossom_set->find(_graph.u(e));
3194 3254
            int right_blossom = _blossom_set->find(_graph.v(e));
3195 3255

	
3196 3256
            if (left_blossom == right_blossom) {
3197 3257
              _delta3->pop();
3198 3258
            } else {
3199 3259
              int left_tree = _tree_set->find(left_blossom);
3200 3260
              int right_tree = _tree_set->find(right_blossom);
3201 3261

	
3202 3262
              if (left_tree == right_tree) {
3203 3263
                shrinkOnEdge(e, left_tree);
3204 3264
              } else {
3205 3265
                augmentOnEdge(e);
3206 3266
                _unmatched -= 2;
3207 3267
              }
3208 3268
            }
3209 3269
          } break;
3210 3270
        case D4:
3211 3271
          splitBlossom(_delta4->top());
3212 3272
          break;
3213 3273
        }
3214 3274
      }
3215 3275
      extractMatching();
3216 3276
      return true;
3217 3277
    }
3218 3278

	
3219 3279
    /// \brief Run the algorithm.
3220 3280
    ///
3221 3281
    /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
3222 3282
    ///
3223 3283
    /// \note mwpm.run() is just a shortcut of the following code.
3224 3284
    /// \code
3225 3285
    ///   mwpm.fractionalInit();
3226 3286
    ///   mwpm.start();
3227 3287
    /// \endcode
3228 3288
    bool run() {
3229 3289
      fractionalInit();
3230 3290
      return start();
3231 3291
    }
3232 3292

	
3233 3293
    /// @}
3234 3294

	
3235 3295
    /// \name Primal Solution
3236 3296
    /// Functions to get the primal solution, i.e. the maximum weighted
3237 3297
    /// perfect matching.\n
3238 3298
    /// Either \ref run() or \ref start() function should be called before
3239 3299
    /// using them.
3240 3300

	
3241 3301
    /// @{
3242 3302

	
3243 3303
    /// \brief Return the weight of the matching.
3244 3304
    ///
3245 3305
    /// This function returns the weight of the found matching.
3246 3306
    ///
3247 3307
    /// \pre Either run() or start() must be called before using this function.
3248 3308
    Value matchingWeight() const {
3249 3309
      Value sum = 0;
3250 3310
      for (NodeIt n(_graph); n != INVALID; ++n) {
3251 3311
        if ((*_matching)[n] != INVALID) {
3252 3312
          sum += _weight[(*_matching)[n]];
3253 3313
        }
3254 3314
      }
3255 3315
      return sum / 2;
3256 3316
    }
3257 3317

	
3258 3318
    /// \brief Return \c true if the given edge is in the matching.
3259 3319
    ///
3260 3320
    /// This function returns \c true if the given edge is in the found
3261 3321
    /// matching.
3262 3322
    ///
3263 3323
    /// \pre Either run() or start() must be called before using this function.
3264 3324
    bool matching(const Edge& edge) const {
3265 3325
      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
3266 3326
    }
3267 3327

	
3268 3328
    /// \brief Return the matching arc (or edge) incident to the given node.
3269 3329
    ///
3270 3330
    /// This function returns the matching arc (or edge) incident to the
3271 3331
    /// given node in the found matching or \c INVALID if the node is
3272 3332
    /// not covered by the matching.
3273 3333
    ///
3274 3334
    /// \pre Either run() or start() must be called before using this function.
3275 3335
    Arc matching(const Node& node) const {
3276 3336
      return (*_matching)[node];
3277 3337
    }
3278 3338

	
3279 3339
    /// \brief Return a const reference to the matching map.
3280 3340
    ///
3281 3341
    /// This function returns a const reference to a node map that stores
3282 3342
    /// the matching arc (or edge) incident to each node.
3283 3343
    const MatchingMap& matchingMap() const {
3284 3344
      return *_matching;
3285 3345
    }
3286 3346

	
3287 3347
    /// \brief Return the mate of the given node.
3288 3348
    ///
3289 3349
    /// This function returns the mate of the given node in the found
3290 3350
    /// matching or \c INVALID if the node is not covered by the matching.
3291 3351
    ///
3292 3352
    /// \pre Either run() or start() must be called before using this function.
3293 3353
    Node mate(const Node& node) const {
3294 3354
      return _graph.target((*_matching)[node]);
3295 3355
    }
3296 3356

	
3297 3357
    /// @}
3298 3358

	
3299 3359
    /// \name Dual Solution
3300 3360
    /// Functions to get the dual solution.\n
3301 3361
    /// Either \ref run() or \ref start() function should be called before
3302 3362
    /// using them.
3303 3363

	
3304 3364
    /// @{
3305 3365

	
3306 3366
    /// \brief Return the value of the dual solution.
3307 3367
    ///
3308 3368
    /// This function returns the value of the dual solution.
3309 3369
    /// It should be equal to the primal value scaled by \ref dualScale
3310 3370
    /// "dual scale".
3311 3371
    ///
3312 3372
    /// \pre Either run() or start() must be called before using this function.
3313 3373
    Value dualValue() const {
3314 3374
      Value sum = 0;
3315 3375
      for (NodeIt n(_graph); n != INVALID; ++n) {
3316 3376
        sum += nodeValue(n);
3317 3377
      }
3318 3378
      for (int i = 0; i < blossomNum(); ++i) {
3319 3379
        sum += blossomValue(i) * (blossomSize(i) / 2);
3320 3380
      }
3321 3381
      return sum;
3322 3382
    }
3323 3383

	
3324 3384
    /// \brief Return the dual value (potential) of the given node.
3325 3385
    ///
3326 3386
    /// This function returns the dual value (potential) of the given node.
3327 3387
    ///
3328 3388
    /// \pre Either run() or start() must be called before using this function.
3329 3389
    Value nodeValue(const Node& n) const {
3330 3390
      return (*_node_potential)[n];
3331 3391
    }
3332 3392

	
3333 3393
    /// \brief Return the number of the blossoms in the basis.
3334 3394
    ///
3335 3395
    /// This function returns the number of the blossoms in the basis.
3336 3396
    ///
3337 3397
    /// \pre Either run() or start() must be called before using this function.
3338 3398
    /// \see BlossomIt
3339 3399
    int blossomNum() const {
3340 3400
      return _blossom_potential.size();
3341 3401
    }
3342 3402

	
3343 3403
    /// \brief Return the number of the nodes in the given blossom.
3344 3404
    ///
3345 3405
    /// This function returns the number of the nodes in the given blossom.
3346 3406
    ///
3347 3407
    /// \pre Either run() or start() must be called before using this function.
3348 3408
    /// \see BlossomIt
3349 3409
    int blossomSize(int k) const {
3350 3410
      return _blossom_potential[k].end - _blossom_potential[k].begin;
3351 3411
    }
3352 3412

	
3353 3413
    /// \brief Return the dual value (ptential) of the given blossom.
3354 3414
    ///
3355 3415
    /// This function returns the dual value (ptential) of the given blossom.
3356 3416
    ///
3357 3417
    /// \pre Either run() or start() must be called before using this function.
3358 3418
    Value blossomValue(int k) const {
3359 3419
      return _blossom_potential[k].value;
3360 3420
    }
3361 3421

	
3362 3422
    /// \brief Iterator for obtaining the nodes of a blossom.
3363 3423
    ///
3364 3424
    /// This class provides an iterator for obtaining the nodes of the
3365 3425
    /// given blossom. It lists a subset of the nodes.
3366 3426
    /// Before using this iterator, you must allocate a
3367 3427
    /// MaxWeightedPerfectMatching class and execute it.
3368 3428
    class BlossomIt {
3369 3429
    public:
3370 3430

	
3371 3431
      /// \brief Constructor.
3372 3432
      ///
3373 3433
      /// Constructor to get the nodes of the given variable.
3374 3434
      ///
3375 3435
      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
3376 3436
      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
3377 3437
      /// must be called before initializing this iterator.
Ignore white space 768 line context
... ...
@@ -907,768 +907,777 @@
907 907
      ///
908 908
      /// It converts the iterator to the current item.
909 909
      operator const Item&() const {
910 910
        return extendFind->items[idx].item;
911 911
      }
912 912

	
913 913
      /// \brief Equality operator
914 914
      ///
915 915
      /// Equality operator
916 916
      bool operator==(const ItemIt& i) {
917 917
        return i.idx == idx;
918 918
      }
919 919

	
920 920
      /// \brief Inequality operator
921 921
      ///
922 922
      /// Inequality operator
923 923
      bool operator!=(const ItemIt& i) {
924 924
        return i.idx != idx;
925 925
      }
926 926

	
927 927
    private:
928 928
      const ExtendFindEnum* extendFind;
929 929
      int idx, fdx;
930 930
    };
931 931

	
932 932
  };
933 933

	
934 934
  /// \ingroup auxdat
935 935
  ///
936 936
  /// \brief A \e Union-Find data structure implementation which
937 937
  /// is able to store a priority for each item and retrieve the minimum of
938 938
  /// each class.
939 939
  ///
940 940
  /// A \e Union-Find data structure implementation which is able to
941 941
  /// store a priority for each item and retrieve the minimum of each
942 942
  /// class. In addition, it supports the joining and splitting the
943 943
  /// components. If you don't need this feature then you makes
944 944
  /// better to use the \ref UnionFind class which is more efficient.
945 945
  ///
946 946
  /// The union-find data strcuture based on a (2, 16)-tree with a
947 947
  /// tournament minimum selection on the internal nodes. The insert
948 948
  /// operation takes O(1), the find, set, decrease and increase takes
949 949
  /// O(log(n)), where n is the number of nodes in the current
950 950
  /// component.  The complexity of join and split is O(log(n)*k),
951 951
  /// where n is the sum of the number of the nodes and k is the
952 952
  /// number of joined components or the number of the components
953 953
  /// after the split.
954 954
  ///
955 955
  /// \pre You need to add all the elements by the \ref insert()
956 956
  /// method.
957 957
  template <typename V, typename IM, typename Comp = std::less<V> >
958 958
  class HeapUnionFind {
959 959
  public:
960 960

	
961 961
    ///\e
962 962
    typedef V Value;
963 963
    ///\e
964 964
    typedef typename IM::Key Item;
965 965
    ///\e
966 966
    typedef IM ItemIntMap;
967 967
    ///\e
968 968
    typedef Comp Compare;
969 969

	
970 970
  private:
971 971

	
972 972
    static const int cmax = 16;
973 973

	
974 974
    ItemIntMap& index;
975 975

	
976 976
    struct ClassNode {
977 977
      int parent;
978 978
      int depth;
979 979

	
980 980
      int left, right;
981 981
      int next, prev;
982 982
    };
983 983

	
984 984
    int first_class;
985 985
    int first_free_class;
986 986
    std::vector<ClassNode> classes;
987 987

	
988 988
    int newClass() {
989 989
      if (first_free_class < 0) {
990 990
        int id = classes.size();
991 991
        classes.push_back(ClassNode());
992 992
        return id;
993 993
      } else {
994 994
        int id = first_free_class;
995 995
        first_free_class = classes[id].next;
996 996
        return id;
997 997
      }
998 998
    }
999 999

	
1000 1000
    void deleteClass(int id) {
1001 1001
      classes[id].next = first_free_class;
1002 1002
      first_free_class = id;
1003 1003
    }
1004 1004

	
1005 1005
    struct ItemNode {
1006 1006
      int parent;
1007 1007
      Item item;
1008 1008
      Value prio;
1009 1009
      int next, prev;
1010 1010
      int left, right;
1011 1011
      int size;
1012 1012
    };
1013 1013

	
1014 1014
    int first_free_node;
1015 1015
    std::vector<ItemNode> nodes;
1016 1016

	
1017 1017
    int newNode() {
1018 1018
      if (first_free_node < 0) {
1019 1019
        int id = nodes.size();
1020 1020
        nodes.push_back(ItemNode());
1021 1021
        return id;
1022 1022
      } else {
1023 1023
        int id = first_free_node;
1024 1024
        first_free_node = nodes[id].next;
1025 1025
        return id;
1026 1026
      }
1027 1027
    }
1028 1028

	
1029 1029
    void deleteNode(int id) {
1030 1030
      nodes[id].next = first_free_node;
1031 1031
      first_free_node = id;
1032 1032
    }
1033 1033

	
1034 1034
    Comp comp;
1035 1035

	
1036 1036
    int findClass(int id) const {
1037 1037
      int kd = id;
1038 1038
      while (kd >= 0) {
1039 1039
        kd = nodes[kd].parent;
1040 1040
      }
1041 1041
      return ~kd;
1042 1042
    }
1043 1043

	
1044 1044
    int leftNode(int id) const {
1045 1045
      int kd = ~(classes[id].parent);
1046 1046
      for (int i = 0; i < classes[id].depth; ++i) {
1047 1047
        kd = nodes[kd].left;
1048 1048
      }
1049 1049
      return kd;
1050 1050
    }
1051 1051

	
1052 1052
    int nextNode(int id) const {
1053 1053
      int depth = 0;
1054 1054
      while (id >= 0 && nodes[id].next == -1) {
1055 1055
        id = nodes[id].parent;
1056 1056
        ++depth;
1057 1057
      }
1058 1058
      if (id < 0) {
1059 1059
        return -1;
1060 1060
      }
1061 1061
      id = nodes[id].next;
1062 1062
      while (depth--) {
1063 1063
        id = nodes[id].left;
1064 1064
      }
1065 1065
      return id;
1066 1066
    }
1067 1067

	
1068 1068

	
1069 1069
    void setPrio(int id) {
1070 1070
      int jd = nodes[id].left;
1071 1071
      nodes[id].prio = nodes[jd].prio;
1072 1072
      nodes[id].item = nodes[jd].item;
1073 1073
      jd = nodes[jd].next;
1074 1074
      while (jd != -1) {
1075 1075
        if (comp(nodes[jd].prio, nodes[id].prio)) {
1076 1076
          nodes[id].prio = nodes[jd].prio;
1077 1077
          nodes[id].item = nodes[jd].item;
1078 1078
        }
1079 1079
        jd = nodes[jd].next;
1080 1080
      }
1081 1081
    }
1082 1082

	
1083 1083
    void push(int id, int jd) {
1084 1084
      nodes[id].size = 1;
1085 1085
      nodes[id].left = nodes[id].right = jd;
1086 1086
      nodes[jd].next = nodes[jd].prev = -1;
1087 1087
      nodes[jd].parent = id;
1088 1088
    }
1089 1089

	
1090 1090
    void pushAfter(int id, int jd) {
1091 1091
      int kd = nodes[id].parent;
1092 1092
      if (nodes[id].next != -1) {
1093 1093
        nodes[nodes[id].next].prev = jd;
1094 1094
        if (kd >= 0) {
1095 1095
          nodes[kd].size += 1;
1096 1096
        }
1097 1097
      } else {
1098 1098
        if (kd >= 0) {
1099 1099
          nodes[kd].right = jd;
1100 1100
          nodes[kd].size += 1;
1101 1101
        }
1102 1102
      }
1103 1103
      nodes[jd].next = nodes[id].next;
1104 1104
      nodes[jd].prev = id;
1105 1105
      nodes[id].next = jd;
1106 1106
      nodes[jd].parent = kd;
1107 1107
    }
1108 1108

	
1109 1109
    void pushRight(int id, int jd) {
1110 1110
      nodes[id].size += 1;
1111 1111
      nodes[jd].prev = nodes[id].right;
1112 1112
      nodes[jd].next = -1;
1113 1113
      nodes[nodes[id].right].next = jd;
1114 1114
      nodes[id].right = jd;
1115 1115
      nodes[jd].parent = id;
1116 1116
    }
1117 1117

	
1118 1118
    void popRight(int id) {
1119 1119
      nodes[id].size -= 1;
1120 1120
      int jd = nodes[id].right;
1121 1121
      nodes[nodes[jd].prev].next = -1;
1122 1122
      nodes[id].right = nodes[jd].prev;
1123 1123
    }
1124 1124

	
1125 1125
    void splice(int id, int jd) {
1126 1126
      nodes[id].size += nodes[jd].size;
1127 1127
      nodes[nodes[id].right].next = nodes[jd].left;
1128 1128
      nodes[nodes[jd].left].prev = nodes[id].right;
1129 1129
      int kd = nodes[jd].left;
1130 1130
      while (kd != -1) {
1131 1131
        nodes[kd].parent = id;
1132 1132
        kd = nodes[kd].next;
1133 1133
      }
1134 1134
      nodes[id].right = nodes[jd].right;
1135 1135
    }
1136 1136

	
1137 1137
    void split(int id, int jd) {
1138 1138
      int kd = nodes[id].parent;
1139 1139
      nodes[kd].right = nodes[id].prev;
1140 1140
      nodes[nodes[id].prev].next = -1;
1141 1141

	
1142 1142
      nodes[jd].left = id;
1143 1143
      nodes[id].prev = -1;
1144 1144
      int num = 0;
1145 1145
      while (id != -1) {
1146 1146
        nodes[id].parent = jd;
1147 1147
        nodes[jd].right = id;
1148 1148
        id = nodes[id].next;
1149 1149
        ++num;
1150 1150
      }
1151 1151
      nodes[kd].size -= num;
1152 1152
      nodes[jd].size = num;
1153 1153
    }
1154 1154

	
1155 1155
    void pushLeft(int id, int jd) {
1156 1156
      nodes[id].size += 1;
1157 1157
      nodes[jd].next = nodes[id].left;
1158 1158
      nodes[jd].prev = -1;
1159 1159
      nodes[nodes[id].left].prev = jd;
1160 1160
      nodes[id].left = jd;
1161 1161
      nodes[jd].parent = id;
1162 1162
    }
1163 1163

	
1164 1164
    void popLeft(int id) {
1165 1165
      nodes[id].size -= 1;
1166 1166
      int jd = nodes[id].left;
1167 1167
      nodes[nodes[jd].next].prev = -1;
1168 1168
      nodes[id].left = nodes[jd].next;
1169 1169
    }
1170 1170

	
1171 1171
    void repairLeft(int id) {
1172 1172
      int jd = ~(classes[id].parent);
1173 1173
      while (nodes[jd].left != -1) {
1174 1174
        int kd = nodes[jd].left;
1175 1175
        if (nodes[jd].size == 1) {
1176 1176
          if (nodes[jd].parent < 0) {
1177 1177
            classes[id].parent = ~kd;
1178 1178
            classes[id].depth -= 1;
1179 1179
            nodes[kd].parent = ~id;
1180 1180
            deleteNode(jd);
1181 1181
            jd = kd;
1182 1182
          } else {
1183 1183
            int pd = nodes[jd].parent;
1184 1184
            if (nodes[nodes[jd].next].size < cmax) {
1185 1185
              pushLeft(nodes[jd].next, nodes[jd].left);
1186 1186
              if (less(jd, nodes[jd].next) ||
1187 1187
                  nodes[jd].item == nodes[pd].item) {
1188 1188
                nodes[nodes[jd].next].prio = nodes[jd].prio;
1189 1189
                nodes[nodes[jd].next].item = nodes[jd].item;
1190 1190
              }
1191 1191
              popLeft(pd);
1192 1192
              deleteNode(jd);
1193 1193
              jd = pd;
1194 1194
            } else {
1195 1195
              int ld = nodes[nodes[jd].next].left;
1196 1196
              popLeft(nodes[jd].next);
1197 1197
              pushRight(jd, ld);
1198 1198
              if (less(ld, nodes[jd].left) ||
1199 1199
                  nodes[ld].item == nodes[pd].item) {
1200 1200
                nodes[jd].item = nodes[ld].item;
1201 1201
                nodes[jd].prio = nodes[ld].prio;
1202 1202
              }
1203 1203
              if (nodes[nodes[jd].next].item == nodes[ld].item) {
1204 1204
                setPrio(nodes[jd].next);
1205 1205
              }
1206 1206
              jd = nodes[jd].left;
1207 1207
            }
1208 1208
          }
1209 1209
        } else {
1210 1210
          jd = nodes[jd].left;
1211 1211
        }
1212 1212
      }
1213 1213
    }
1214 1214

	
1215 1215
    void repairRight(int id) {
1216 1216
      int jd = ~(classes[id].parent);
1217 1217
      while (nodes[jd].right != -1) {
1218 1218
        int kd = nodes[jd].right;
1219 1219
        if (nodes[jd].size == 1) {
1220 1220
          if (nodes[jd].parent < 0) {
1221 1221
            classes[id].parent = ~kd;
1222 1222
            classes[id].depth -= 1;
1223 1223
            nodes[kd].parent = ~id;
1224 1224
            deleteNode(jd);
1225 1225
            jd = kd;
1226 1226
          } else {
1227 1227
            int pd = nodes[jd].parent;
1228 1228
            if (nodes[nodes[jd].prev].size < cmax) {
1229 1229
              pushRight(nodes[jd].prev, nodes[jd].right);
1230 1230
              if (less(jd, nodes[jd].prev) ||
1231 1231
                  nodes[jd].item == nodes[pd].item) {
1232 1232
                nodes[nodes[jd].prev].prio = nodes[jd].prio;
1233 1233
                nodes[nodes[jd].prev].item = nodes[jd].item;
1234 1234
              }
1235 1235
              popRight(pd);
1236 1236
              deleteNode(jd);
1237 1237
              jd = pd;
1238 1238
            } else {
1239 1239
              int ld = nodes[nodes[jd].prev].right;
1240 1240
              popRight(nodes[jd].prev);
1241 1241
              pushLeft(jd, ld);
1242 1242
              if (less(ld, nodes[jd].right) ||
1243 1243
                  nodes[ld].item == nodes[pd].item) {
1244 1244
                nodes[jd].item = nodes[ld].item;
1245 1245
                nodes[jd].prio = nodes[ld].prio;
1246 1246
              }
1247 1247
              if (nodes[nodes[jd].prev].item == nodes[ld].item) {
1248 1248
                setPrio(nodes[jd].prev);
1249 1249
              }
1250 1250
              jd = nodes[jd].right;
1251 1251
            }
1252 1252
          }
1253 1253
        } else {
1254 1254
          jd = nodes[jd].right;
1255 1255
        }
1256 1256
      }
1257 1257
    }
1258 1258

	
1259 1259

	
1260 1260
    bool less(int id, int jd) const {
1261 1261
      return comp(nodes[id].prio, nodes[jd].prio);
1262 1262
    }
1263 1263

	
1264 1264
  public:
1265 1265

	
1266 1266
    /// \brief Returns true when the given class is alive.
1267 1267
    ///
1268 1268
    /// Returns true when the given class is alive, ie. the class is
1269 1269
    /// not nested into other class.
1270 1270
    bool alive(int cls) const {
1271 1271
      return classes[cls].parent < 0;
1272 1272
    }
1273 1273

	
1274 1274
    /// \brief Returns true when the given class is trivial.
1275 1275
    ///
1276 1276
    /// Returns true when the given class is trivial, ie. the class
1277 1277
    /// contains just one item directly.
1278 1278
    bool trivial(int cls) const {
1279 1279
      return classes[cls].left == -1;
1280 1280
    }
1281 1281

	
1282 1282
    /// \brief Constructs the union-find.
1283 1283
    ///
1284 1284
    /// Constructs the union-find.
1285 1285
    /// \brief _index The index map of the union-find. The data
1286 1286
    /// structure uses internally for store references.
1287 1287
    HeapUnionFind(ItemIntMap& _index)
1288 1288
      : index(_index), first_class(-1),
1289 1289
        first_free_class(-1), first_free_node(-1) {}
1290 1290

	
1291
    /// \brief Clears the union-find data structure
1292
    ///
1293
    /// Erase each item from the data structure.
1294
    void clear() {
1295
      nodes.clear();
1296
      classes.clear();
1297
      first_free_node = first_free_class = first_class = -1;
1298
    }
1299

	
1291 1300
    /// \brief Insert a new node into a new component.
1292 1301
    ///
1293 1302
    /// Insert a new node into a new component.
1294 1303
    /// \param item The item of the new node.
1295 1304
    /// \param prio The priority of the new node.
1296 1305
    /// \return The class id of the one-item-heap.
1297 1306
    int insert(const Item& item, const Value& prio) {
1298 1307
      int id = newNode();
1299 1308
      nodes[id].item = item;
1300 1309
      nodes[id].prio = prio;
1301 1310
      nodes[id].size = 0;
1302 1311

	
1303 1312
      nodes[id].prev = -1;
1304 1313
      nodes[id].next = -1;
1305 1314

	
1306 1315
      nodes[id].left = -1;
1307 1316
      nodes[id].right = -1;
1308 1317

	
1309 1318
      nodes[id].item = item;
1310 1319
      index[item] = id;
1311 1320

	
1312 1321
      int class_id = newClass();
1313 1322
      classes[class_id].parent = ~id;
1314 1323
      classes[class_id].depth = 0;
1315 1324

	
1316 1325
      classes[class_id].left = -1;
1317 1326
      classes[class_id].right = -1;
1318 1327

	
1319 1328
      if (first_class != -1) {
1320 1329
        classes[first_class].prev = class_id;
1321 1330
      }
1322 1331
      classes[class_id].next = first_class;
1323 1332
      classes[class_id].prev = -1;
1324 1333
      first_class = class_id;
1325 1334

	
1326 1335
      nodes[id].parent = ~class_id;
1327 1336

	
1328 1337
      return class_id;
1329 1338
    }
1330 1339

	
1331 1340
    /// \brief The class of the item.
1332 1341
    ///
1333 1342
    /// \return The alive class id of the item, which is not nested into
1334 1343
    /// other classes.
1335 1344
    ///
1336 1345
    /// The time complexity is O(log(n)).
1337 1346
    int find(const Item& item) const {
1338 1347
      return findClass(index[item]);
1339 1348
    }
1340 1349

	
1341 1350
    /// \brief Joins the classes.
1342 1351
    ///
1343 1352
    /// The current function joins the given classes. The parameter is
1344 1353
    /// an STL range which should be contains valid class ids. The
1345 1354
    /// time complexity is O(log(n)*k) where n is the overall number
1346 1355
    /// of the joined nodes and k is the number of classes.
1347 1356
    /// \return The class of the joined classes.
1348 1357
    /// \pre The range should contain at least two class ids.
1349 1358
    template <typename Iterator>
1350 1359
    int join(Iterator begin, Iterator end) {
1351 1360
      std::vector<int> cs;
1352 1361
      for (Iterator it = begin; it != end; ++it) {
1353 1362
        cs.push_back(*it);
1354 1363
      }
1355 1364

	
1356 1365
      int class_id = newClass();
1357 1366
      { // creation union-find
1358 1367

	
1359 1368
        if (first_class != -1) {
1360 1369
          classes[first_class].prev = class_id;
1361 1370
        }
1362 1371
        classes[class_id].next = first_class;
1363 1372
        classes[class_id].prev = -1;
1364 1373
        first_class = class_id;
1365 1374

	
1366 1375
        classes[class_id].depth = classes[cs[0]].depth;
1367 1376
        classes[class_id].parent = classes[cs[0]].parent;
1368 1377
        nodes[~(classes[class_id].parent)].parent = ~class_id;
1369 1378

	
1370 1379
        int l = cs[0];
1371 1380

	
1372 1381
        classes[class_id].left = l;
1373 1382
        classes[class_id].right = l;
1374 1383

	
1375 1384
        if (classes[l].next != -1) {
1376 1385
          classes[classes[l].next].prev = classes[l].prev;
1377 1386
        }
1378 1387
        classes[classes[l].prev].next = classes[l].next;
1379 1388

	
1380 1389
        classes[l].prev = -1;
1381 1390
        classes[l].next = -1;
1382 1391

	
1383 1392
        classes[l].depth = leftNode(l);
1384 1393
        classes[l].parent = class_id;
1385 1394

	
1386 1395
      }
1387 1396

	
1388 1397
      { // merging of heap
1389 1398
        int l = class_id;
1390 1399
        for (int ci = 1; ci < int(cs.size()); ++ci) {
1391 1400
          int r = cs[ci];
1392 1401
          int rln = leftNode(r);
1393 1402
          if (classes[l].depth > classes[r].depth) {
1394 1403
            int id = ~(classes[l].parent);
1395 1404
            for (int i = classes[r].depth + 1; i < classes[l].depth; ++i) {
1396 1405
              id = nodes[id].right;
1397 1406
            }
1398 1407
            while (id >= 0 && nodes[id].size == cmax) {
1399 1408
              int new_id = newNode();
1400 1409
              int right_id = nodes[id].right;
1401 1410

	
1402 1411
              popRight(id);
1403 1412
              if (nodes[id].item == nodes[right_id].item) {
1404 1413
                setPrio(id);
1405 1414
              }
1406 1415
              push(new_id, right_id);
1407 1416
              pushRight(new_id, ~(classes[r].parent));
1408 1417

	
1409 1418
              if (less(~classes[r].parent, right_id)) {
1410 1419
                nodes[new_id].item = nodes[~classes[r].parent].item;
1411 1420
                nodes[new_id].prio = nodes[~classes[r].parent].prio;
1412 1421
              } else {
1413 1422
                nodes[new_id].item = nodes[right_id].item;
1414 1423
                nodes[new_id].prio = nodes[right_id].prio;
1415 1424
              }
1416 1425

	
1417 1426
              id = nodes[id].parent;
1418 1427
              classes[r].parent = ~new_id;
1419 1428
            }
1420 1429
            if (id < 0) {
1421 1430
              int new_parent = newNode();
1422 1431
              nodes[new_parent].next = -1;
1423 1432
              nodes[new_parent].prev = -1;
1424 1433
              nodes[new_parent].parent = ~l;
1425 1434

	
1426 1435
              push(new_parent, ~(classes[l].parent));
1427 1436
              pushRight(new_parent, ~(classes[r].parent));
1428 1437
              setPrio(new_parent);
1429 1438

	
1430 1439
              classes[l].parent = ~new_parent;
1431 1440
              classes[l].depth += 1;
1432 1441
            } else {
1433 1442
              pushRight(id, ~(classes[r].parent));
1434 1443
              while (id >= 0 && less(~(classes[r].parent), id)) {
1435 1444
                nodes[id].prio = nodes[~(classes[r].parent)].prio;
1436 1445
                nodes[id].item = nodes[~(classes[r].parent)].item;
1437 1446
                id = nodes[id].parent;
1438 1447
              }
1439 1448
            }
1440 1449
          } else if (classes[r].depth > classes[l].depth) {
1441 1450
            int id = ~(classes[r].parent);
1442 1451
            for (int i = classes[l].depth + 1; i < classes[r].depth; ++i) {
1443 1452
              id = nodes[id].left;
1444 1453
            }
1445 1454
            while (id >= 0 && nodes[id].size == cmax) {
1446 1455
              int new_id = newNode();
1447 1456
              int left_id = nodes[id].left;
1448 1457

	
1449 1458
              popLeft(id);
1450 1459
              if (nodes[id].prio == nodes[left_id].prio) {
1451 1460
                setPrio(id);
1452 1461
              }
1453 1462
              push(new_id, left_id);
1454 1463
              pushLeft(new_id, ~(classes[l].parent));
1455 1464

	
1456 1465
              if (less(~classes[l].parent, left_id)) {
1457 1466
                nodes[new_id].item = nodes[~classes[l].parent].item;
1458 1467
                nodes[new_id].prio = nodes[~classes[l].parent].prio;
1459 1468
              } else {
1460 1469
                nodes[new_id].item = nodes[left_id].item;
1461 1470
                nodes[new_id].prio = nodes[left_id].prio;
1462 1471
              }
1463 1472

	
1464 1473
              id = nodes[id].parent;
1465 1474
              classes[l].parent = ~new_id;
1466 1475

	
1467 1476
            }
1468 1477
            if (id < 0) {
1469 1478
              int new_parent = newNode();
1470 1479
              nodes[new_parent].next = -1;
1471 1480
              nodes[new_parent].prev = -1;
1472 1481
              nodes[new_parent].parent = ~l;
1473 1482

	
1474 1483
              push(new_parent, ~(classes[r].parent));
1475 1484
              pushLeft(new_parent, ~(classes[l].parent));
1476 1485
              setPrio(new_parent);
1477 1486

	
1478 1487
              classes[r].parent = ~new_parent;
1479 1488
              classes[r].depth += 1;
1480 1489
            } else {
1481 1490
              pushLeft(id, ~(classes[l].parent));
1482 1491
              while (id >= 0 && less(~(classes[l].parent), id)) {
1483 1492
                nodes[id].prio = nodes[~(classes[l].parent)].prio;
1484 1493
                nodes[id].item = nodes[~(classes[l].parent)].item;
1485 1494
                id = nodes[id].parent;
1486 1495
              }
1487 1496
            }
1488 1497
            nodes[~(classes[r].parent)].parent = ~l;
1489 1498
            classes[l].parent = classes[r].parent;
1490 1499
            classes[l].depth = classes[r].depth;
1491 1500
          } else {
1492 1501
            if (classes[l].depth != 0 &&
1493 1502
                nodes[~(classes[l].parent)].size +
1494 1503
                nodes[~(classes[r].parent)].size <= cmax) {
1495 1504
              splice(~(classes[l].parent), ~(classes[r].parent));
1496 1505
              deleteNode(~(classes[r].parent));
1497 1506
              if (less(~(classes[r].parent), ~(classes[l].parent))) {
1498 1507
                nodes[~(classes[l].parent)].prio =
1499 1508
                  nodes[~(classes[r].parent)].prio;
1500 1509
                nodes[~(classes[l].parent)].item =
1501 1510
                  nodes[~(classes[r].parent)].item;
1502 1511
              }
1503 1512
            } else {
1504 1513
              int new_parent = newNode();
1505 1514
              nodes[new_parent].next = nodes[new_parent].prev = -1;
1506 1515
              push(new_parent, ~(classes[l].parent));
1507 1516
              pushRight(new_parent, ~(classes[r].parent));
1508 1517
              setPrio(new_parent);
1509 1518

	
1510 1519
              classes[l].parent = ~new_parent;
1511 1520
              classes[l].depth += 1;
1512 1521
              nodes[new_parent].parent = ~l;
1513 1522
            }
1514 1523
          }
1515 1524
          if (classes[r].next != -1) {
1516 1525
            classes[classes[r].next].prev = classes[r].prev;
1517 1526
          }
1518 1527
          classes[classes[r].prev].next = classes[r].next;
1519 1528

	
1520 1529
          classes[r].prev = classes[l].right;
1521 1530
          classes[classes[l].right].next = r;
1522 1531
          classes[l].right = r;
1523 1532
          classes[r].parent = l;
1524 1533

	
1525 1534
          classes[r].next = -1;
1526 1535
          classes[r].depth = rln;
1527 1536
        }
1528 1537
      }
1529 1538
      return class_id;
1530 1539
    }
1531 1540

	
1532 1541
    /// \brief Split the class to subclasses.
1533 1542
    ///
1534 1543
    /// The current function splits the given class. The join, which
1535 1544
    /// made the current class, stored a reference to the
1536 1545
    /// subclasses. The \c splitClass() member restores the classes
1537 1546
    /// and creates the heaps. The parameter is an STL output iterator
1538 1547
    /// which will be filled with the subclass ids. The time
1539 1548
    /// complexity is O(log(n)*k) where n is the overall number of
1540 1549
    /// nodes in the splitted classes and k is the number of the
1541 1550
    /// classes.
1542 1551
    template <typename Iterator>
1543 1552
    void split(int cls, Iterator out) {
1544 1553
      std::vector<int> cs;
1545 1554
      { // splitting union-find
1546 1555
        int id = cls;
1547 1556
        int l = classes[id].left;
1548 1557

	
1549 1558
        classes[l].parent = classes[id].parent;
1550 1559
        classes[l].depth = classes[id].depth;
1551 1560

	
1552 1561
        nodes[~(classes[l].parent)].parent = ~l;
1553 1562

	
1554 1563
        *out++ = l;
1555 1564

	
1556 1565
        while (l != -1) {
1557 1566
          cs.push_back(l);
1558 1567
          l = classes[l].next;
1559 1568
        }
1560 1569

	
1561 1570
        classes[classes[id].right].next = first_class;
1562 1571
        classes[first_class].prev = classes[id].right;
1563 1572
        first_class = classes[id].left;
1564 1573

	
1565 1574
        if (classes[id].next != -1) {
1566 1575
          classes[classes[id].next].prev = classes[id].prev;
1567 1576
        }
1568 1577
        classes[classes[id].prev].next = classes[id].next;
1569 1578

	
1570 1579
        deleteClass(id);
1571 1580
      }
1572 1581

	
1573 1582
      {
1574 1583
        for (int i = 1; i < int(cs.size()); ++i) {
1575 1584
          int l = classes[cs[i]].depth;
1576 1585
          while (nodes[nodes[l].parent].left == l) {
1577 1586
            l = nodes[l].parent;
1578 1587
          }
1579 1588
          int r = l;
1580 1589
          while (nodes[l].parent >= 0) {
1581 1590
            l = nodes[l].parent;
1582 1591
            int new_node = newNode();
1583 1592

	
1584 1593
            nodes[new_node].prev = -1;
1585 1594
            nodes[new_node].next = -1;
1586 1595

	
1587 1596
            split(r, new_node);
1588 1597
            pushAfter(l, new_node);
1589 1598
            setPrio(l);
1590 1599
            setPrio(new_node);
1591 1600
            r = new_node;
1592 1601
          }
1593 1602
          classes[cs[i]].parent = ~r;
1594 1603
          classes[cs[i]].depth = classes[~(nodes[l].parent)].depth;
1595 1604
          nodes[r].parent = ~cs[i];
1596 1605

	
1597 1606
          nodes[l].next = -1;
1598 1607
          nodes[r].prev = -1;
1599 1608

	
1600 1609
          repairRight(~(nodes[l].parent));
1601 1610
          repairLeft(cs[i]);
1602 1611

	
1603 1612
          *out++ = cs[i];
1604 1613
        }
1605 1614
      }
1606 1615
    }
1607 1616

	
1608 1617
    /// \brief Gives back the priority of the current item.
1609 1618
    ///
1610 1619
    /// Gives back the priority of the current item.
1611 1620
    const Value& operator[](const Item& item) const {
1612 1621
      return nodes[index[item]].prio;
1613 1622
    }
1614 1623

	
1615 1624
    /// \brief Sets the priority of the current item.
1616 1625
    ///
1617 1626
    /// Sets the priority of the current item.
1618 1627
    void set(const Item& item, const Value& prio) {
1619 1628
      if (comp(prio, nodes[index[item]].prio)) {
1620 1629
        decrease(item, prio);
1621 1630
      } else if (!comp(prio, nodes[index[item]].prio)) {
1622 1631
        increase(item, prio);
1623 1632
      }
1624 1633
    }
1625 1634

	
1626 1635
    /// \brief Increase the priority of the current item.
1627 1636
    ///
1628 1637
    /// Increase the priority of the current item.
1629 1638
    void increase(const Item& item, const Value& prio) {
1630 1639
      int id = index[item];
1631 1640
      int kd = nodes[id].parent;
1632 1641
      nodes[id].prio = prio;
1633 1642
      while (kd >= 0 && nodes[kd].item == item) {
1634 1643
        setPrio(kd);
1635 1644
        kd = nodes[kd].parent;
1636 1645
      }
1637 1646
    }
1638 1647

	
1639 1648
    /// \brief Increase the priority of the current item.
1640 1649
    ///
1641 1650
    /// Increase the priority of the current item.
1642 1651
    void decrease(const Item& item, const Value& prio) {
1643 1652
      int id = index[item];
1644 1653
      int kd = nodes[id].parent;
1645 1654
      nodes[id].prio = prio;
1646 1655
      while (kd >= 0 && less(id, kd)) {
1647 1656
        nodes[kd].prio = prio;
1648 1657
        nodes[kd].item = item;
1649 1658
        kd = nodes[kd].parent;
1650 1659
      }
1651 1660
    }
1652 1661

	
1653 1662
    /// \brief Gives back the minimum priority of the class.
1654 1663
    ///
1655 1664
    /// Gives back the minimum priority of the class.
1656 1665
    const Value& classPrio(int cls) const {
1657 1666
      return nodes[~(classes[cls].parent)].prio;
1658 1667
    }
1659 1668

	
1660 1669
    /// \brief Gives back the minimum priority item of the class.
1661 1670
    ///
1662 1671
    /// \return Gives back the minimum priority item of the class.
1663 1672
    const Item& classTop(int cls) const {
1664 1673
      return nodes[~(classes[cls].parent)].item;
1665 1674
    }
1666 1675

	
1667 1676
    /// \brief Gives back a representant item of the class.
1668 1677
    ///
1669 1678
    /// Gives back a representant item of the class.
1670 1679
    /// The representant is indpendent from the priorities of the
1671 1680
    /// items.
1672 1681
    const Item& classRep(int id) const {
1673 1682
      int parent = classes[id].parent;
1674 1683
      return nodes[parent >= 0 ? classes[id].depth : leftNode(id)].item;
0 comments (0 inline)