Changeset 2548:a3ba22ebccc6 in lemon0.x
 Timestamp:
 12/28/07 12:00:51 (13 years ago)
 Branch:
 default
 Phase:
 public
 Convert:
 svn:c9d7d8f590d60310b91f818b3a526b0e/lemon/trunk@3425
 Files:

 4 edited
Legend:
 Unmodified
 Added
 Removed

doc/groups.dox
r2530 r2548 319 319 for find matchings in graphs and bipartite graphs. 320 320 321 This group provides some algorithm objects and function 322 to calculate matchings in graphs and bipartite graphs. 321 This group provides some algorithm objects and function to calculate 322 matchings in graphs and bipartite graphs. The general matching problem is 323 finding a subset of the edges which does not shares common endpoints. 324 325 There are several different algorithms for calculate matchings in 326 graphs. The matching problems in bipartite graphs are generally 327 easier than in general graphs. The goal of the matching optimization 328 can be the finding maximum cardinality, maximum weight or minimum cost 329 matching. The search can be constrained to find perfect or 330 maximum cardinality matching. 331 332 Lemon contains the next algorithms: 333  \ref lemon::MaxBipartiteMatching "MaxBipartiteMatching" HopcroftKarp 334 augmenting path algorithm for calculate maximum cardinality matching in 335 bipartite graphs 336  \ref lemon::PrBipartiteMatching "PrBipartiteMatching" PushRelabel 337 algorithm for calculate maximum cardinality matching in bipartite graphs 338  \ref lemon::MaxWeightedBipartiteMatching "MaxWeightedBipartiteMatching" 339 Successive shortest path algorithm for calculate maximum weighted matching 340 and maximum weighted bipartite matching in bipartite graph 341  \ref lemon::MinCostMaxBipartiteMatching "MinCostMaxBipartiteMatching" 342 Successive shortest path algorithm for calculate minimum cost maximum 343 matching in bipartite graph 344  \ref lemon::MaxMatching "MaxMatching" Edmond's blossom shrinking algorithm 345 for calculate maximum cardinality matching in general graph 346  \ref lemon::MaxWeightedMatching "MaxWeightedMatching" Edmond's blossom 347 shrinking algorithm for calculate maximum weighted matching in general 348 graph 349  \ref lemon::MaxWeightedPerfectMatching "MaxWeightedPerfectMatching" 350 Edmond's blossom shrinking algorithm for calculate maximum weighted 351 perfect matching in general graph 323 352 324 353 \image html bipartite_matching.png 
lemon/bin_heap.h
r2547 r2548 53 53 54 54 public: 55 ///\e 55 56 typedef _ItemIntMap ItemIntMap; 57 ///\e 56 58 typedef _Prio Prio; 59 ///\e 57 60 typedef typename ItemIntMap::Key Item; 61 ///\e 58 62 typedef std::pair<Item,Prio> Pair; 63 ///\e 59 64 typedef _Compare Compare; 60 65 … … 322 327 } 323 328 329 /// \brief Replaces an item in the heap. 330 /// 331 /// The \c i item is replaced with \c j item. The \c i item should 332 /// be in the heap, while the \c j should be out of the heap. The 333 /// \c i item will out of the heap and \c j will be in the heap 334 /// with the same prioriority as prevoiusly the \c i item. 335 void replace(const Item& i, const Item& j) { 336 int idx = iim[i]; 337 iim.set(i, iim[j]); 338 iim.set(j, idx); 339 data[idx].first = j; 340 } 341 324 342 }; // class BinHeap 325 343 
lemon/max_matching.h
r2505 r2548 20 20 #define LEMON_MAX_MATCHING_H 21 21 22 #include <vector> 22 23 #include <queue> 24 #include <set> 25 23 26 #include <lemon/bits/invalid.h> 24 27 #include <lemon/unionfind.h> 25 28 #include <lemon/graph_utils.h> 29 #include <lemon/bin_heap.h> 26 30 27 31 ///\ingroup matching … … 32 36 33 37 ///\ingroup matching 34 38 /// 35 39 ///\brief Edmonds' alternating forest maximum matching algorithm. 36 40 /// … … 619 623 620 624 }; 625 626 /// \ingroup matching 627 /// 628 /// \brief Weighted matching in general undirected graphs 629 /// 630 /// This class provides an efficient implementation of Edmond's 631 /// maximum weighted matching algorithm. The implementation is based 632 /// on extensive use of priority queues and provides 633 /// \f$O(nm\log(n))\f$ time complexity. 634 /// 635 /// The maximum weighted matching problem is to find undirected 636 /// edges in the graph with maximum overall weight and no two of 637 /// them shares their endpoints. The problem can be formulated with 638 /// the next linear program: 639 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f] 640 ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert  1}{2} \quad \forall B\in\mathcal{O}\f] 641 /// \f[x_e \ge 0\quad \forall e\in E\f] 642 /// \f[\max \sum_{e\in E}x_ew_e\f] 643 /// where \f$\delta(X)\f$ is the set of edges incident to a node in 644 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in 645 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of 646 /// the nodes. 647 /// 648 /// The algorithm calculates an optimal matching and a proof of the 649 /// optimality. The solution of the dual problem can be used to check 650 /// the result of the algorithm. The dual linear problem is the next: 651 /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f] 652 /// \f[y_u \ge 0 \quad \forall u \in V\f] 653 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] 654 /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert  1}{2}z_B\f] 655 /// 656 /// The algorithm can be executed with \c run() or the \c init() and 657 /// then the \c start() member functions. After it the matching can 658 /// be asked with \c matching() or mate() functions. The dual 659 /// solution can be get with \c nodeValue(), \c blossomNum() and \c 660 /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt 661 /// "BlossomIt" nested class which is able to iterate on the nodes 662 /// of a blossom. If the value type is integral then the dual 663 /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4". 664 /// 665 /// \author Balazs Dezso 666 template <typename _UGraph, 667 typename _WeightMap = typename _UGraph::template UEdgeMap<int> > 668 class MaxWeightedMatching { 669 public: 670 671 typedef _UGraph UGraph; 672 typedef _WeightMap WeightMap; 673 typedef typename WeightMap::Value Value; 674 675 /// \brief Scaling factor for dual solution 676 /// 677 /// Scaling factor for dual solution, it is equal to 4 or 1 678 /// according to the value type. 679 static const int dualScale = 680 std::numeric_limits<Value>::is_integer ? 4 : 1; 681 682 typedef typename UGraph::template NodeMap<typename UGraph::Edge> 683 MatchingMap; 684 685 private: 686 687 UGRAPH_TYPEDEFS(typename UGraph); 688 689 typedef typename UGraph::template NodeMap<Value> NodePotential; 690 typedef std::vector<Node> BlossomNodeList; 691 692 struct BlossomVariable { 693 int begin, end; 694 Value value; 695 696 BlossomVariable(int _begin, int _end, Value _value) 697 : begin(_begin), end(_end), value(_value) {} 698 699 }; 700 701 typedef std::vector<BlossomVariable> BlossomPotential; 702 703 const UGraph& _ugraph; 704 const WeightMap& _weight; 705 706 MatchingMap* _matching; 707 708 NodePotential* _node_potential; 709 710 BlossomPotential _blossom_potential; 711 BlossomNodeList _blossom_node_list; 712 713 int _node_num; 714 int _blossom_num; 715 716 typedef typename UGraph::template NodeMap<int> NodeIntMap; 717 typedef typename UGraph::template EdgeMap<int> EdgeIntMap; 718 typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap; 719 typedef IntegerMap<int> IntIntMap; 720 721 enum Status { 722 EVEN = 1, MATCHED = 0, ODD = 1, UNMATCHED = 2 723 }; 724 725 typedef HeapUnionFind<Value, NodeIntMap> BlossomSet; 726 struct BlossomData { 727 int tree; 728 Status status; 729 Edge pred, next; 730 Value pot, offset; 731 Node base; 732 }; 733 734 NodeIntMap *_blossom_index; 735 BlossomSet *_blossom_set; 736 IntegerMap<BlossomData>* _blossom_data; 737 738 NodeIntMap *_node_index; 739 EdgeIntMap *_node_heap_index; 740 741 struct NodeData { 742 743 NodeData(EdgeIntMap& node_heap_index) 744 : heap(node_heap_index) {} 745 746 int blossom; 747 Value pot; 748 BinHeap<Value, EdgeIntMap> heap; 749 std::map<int, Edge> heap_index; 750 751 int tree; 752 }; 753 754 IntegerMap<NodeData>* _node_data; 755 756 typedef ExtendFindEnum<IntIntMap> TreeSet; 757 758 IntIntMap *_tree_set_index; 759 TreeSet *_tree_set; 760 761 NodeIntMap *_delta1_index; 762 BinHeap<Value, NodeIntMap> *_delta1; 763 764 IntIntMap *_delta2_index; 765 BinHeap<Value, IntIntMap> *_delta2; 766 767 UEdgeIntMap *_delta3_index; 768 BinHeap<Value, UEdgeIntMap> *_delta3; 769 770 IntIntMap *_delta4_index; 771 BinHeap<Value, IntIntMap> *_delta4; 772 773 Value _delta_sum; 774 775 void createStructures() { 776 _node_num = countNodes(_ugraph); 777 _blossom_num = _node_num * 3 / 2; 778 779 if (!_matching) { 780 _matching = new MatchingMap(_ugraph); 781 } 782 if (!_node_potential) { 783 _node_potential = new NodePotential(_ugraph); 784 } 785 if (!_blossom_set) { 786 _blossom_index = new NodeIntMap(_ugraph); 787 _blossom_set = new BlossomSet(*_blossom_index); 788 _blossom_data = new IntegerMap<BlossomData>(_blossom_num); 789 } 790 791 if (!_node_index) { 792 _node_index = new NodeIntMap(_ugraph); 793 _node_heap_index = new EdgeIntMap(_ugraph); 794 _node_data = new IntegerMap<NodeData>(_node_num, 795 NodeData(*_node_heap_index)); 796 } 797 798 if (!_tree_set) { 799 _tree_set_index = new IntIntMap(_blossom_num); 800 _tree_set = new TreeSet(*_tree_set_index); 801 } 802 if (!_delta1) { 803 _delta1_index = new NodeIntMap(_ugraph); 804 _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index); 805 } 806 if (!_delta2) { 807 _delta2_index = new IntIntMap(_blossom_num); 808 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index); 809 } 810 if (!_delta3) { 811 _delta3_index = new UEdgeIntMap(_ugraph); 812 _delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index); 813 } 814 if (!_delta4) { 815 _delta4_index = new IntIntMap(_blossom_num); 816 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index); 817 } 818 } 819 820 void destroyStructures() { 821 _node_num = countNodes(_ugraph); 822 _blossom_num = _node_num * 3 / 2; 823 824 if (_matching) { 825 delete _matching; 826 } 827 if (_node_potential) { 828 delete _node_potential; 829 } 830 if (_blossom_set) { 831 delete _blossom_index; 832 delete _blossom_set; 833 delete _blossom_data; 834 } 835 836 if (_node_index) { 837 delete _node_index; 838 delete _node_heap_index; 839 delete _node_data; 840 } 841 842 if (_tree_set) { 843 delete _tree_set_index; 844 delete _tree_set; 845 } 846 if (_delta1) { 847 delete _delta1_index; 848 delete _delta1; 849 } 850 if (_delta2) { 851 delete _delta2_index; 852 delete _delta2; 853 } 854 if (_delta3) { 855 delete _delta3_index; 856 delete _delta3; 857 } 858 if (_delta4) { 859 delete _delta4_index; 860 delete _delta4; 861 } 862 } 863 864 void matchedToEven(int blossom, int tree) { 865 if (_delta2>state(blossom) == _delta2>IN_HEAP) { 866 _delta2>erase(blossom); 867 } 868 869 if (!_blossom_set>trivial(blossom)) { 870 (*_blossom_data)[blossom].pot = 871 2 * (_delta_sum  (*_blossom_data)[blossom].offset); 872 } 873 874 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 875 n != INVALID; ++n) { 876 877 _blossom_set>increase(n, std::numeric_limits<Value>::max()); 878 int ni = (*_node_index)[n]; 879 880 (*_node_data)[ni].heap.clear(); 881 (*_node_data)[ni].heap_index.clear(); 882 883 (*_node_data)[ni].pot += _delta_sum  (*_blossom_data)[blossom].offset; 884 885 _delta1>push(n, (*_node_data)[ni].pot); 886 887 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { 888 Node v = _ugraph.source(e); 889 int vb = _blossom_set>find(v); 890 int vi = (*_node_index)[v]; 891 892 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot  893 dualScale * _weight[e]; 894 895 if ((*_blossom_data)[vb].status == EVEN) { 896 if (_delta3>state(e) != _delta3>IN_HEAP && blossom != vb) { 897 _delta3>push(e, rw / 2); 898 } 899 } else if ((*_blossom_data)[vb].status == UNMATCHED) { 900 if (_delta3>state(e) != _delta3>IN_HEAP) { 901 _delta3>push(e, rw); 902 } 903 } else { 904 typename std::map<int, Edge>::iterator it = 905 (*_node_data)[vi].heap_index.find(tree); 906 907 if (it != (*_node_data)[vi].heap_index.end()) { 908 if ((*_node_data)[vi].heap[it>second] > rw) { 909 (*_node_data)[vi].heap.replace(it>second, e); 910 (*_node_data)[vi].heap.decrease(e, rw); 911 it>second = e; 912 } 913 } else { 914 (*_node_data)[vi].heap.push(e, rw); 915 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); 916 } 917 918 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { 919 _blossom_set>decrease(v, (*_node_data)[vi].heap.prio()); 920 921 if ((*_blossom_data)[vb].status == MATCHED) { 922 if (_delta2>state(vb) != _delta2>IN_HEAP) { 923 _delta2>push(vb, _blossom_set>classPrio(vb)  924 (*_blossom_data)[vb].offset); 925 } else if ((*_delta2)[vb] > _blossom_set>classPrio(vb)  926 (*_blossom_data)[vb].offset){ 927 _delta2>decrease(vb, _blossom_set>classPrio(vb)  928 (*_blossom_data)[vb].offset); 929 } 930 } 931 } 932 } 933 } 934 } 935 (*_blossom_data)[blossom].offset = 0; 936 } 937 938 void matchedToOdd(int blossom) { 939 if (_delta2>state(blossom) == _delta2>IN_HEAP) { 940 _delta2>erase(blossom); 941 } 942 (*_blossom_data)[blossom].offset += _delta_sum; 943 if (!_blossom_set>trivial(blossom)) { 944 _delta4>push(blossom, (*_blossom_data)[blossom].pot / 2 + 945 (*_blossom_data)[blossom].offset); 946 } 947 } 948 949 void evenToMatched(int blossom, int tree) { 950 if (!_blossom_set>trivial(blossom)) { 951 (*_blossom_data)[blossom].pot += 2 * _delta_sum; 952 } 953 954 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 955 n != INVALID; ++n) { 956 int ni = (*_node_index)[n]; 957 (*_node_data)[ni].pot = _delta_sum; 958 959 _delta1>erase(n); 960 961 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { 962 Node v = _ugraph.source(e); 963 int vb = _blossom_set>find(v); 964 int vi = (*_node_index)[v]; 965 966 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot  967 dualScale * _weight[e]; 968 969 if (vb == blossom) { 970 if (_delta3>state(e) == _delta3>IN_HEAP) { 971 _delta3>erase(e); 972 } 973 } else if ((*_blossom_data)[vb].status == EVEN) { 974 975 if (_delta3>state(e) == _delta3>IN_HEAP) { 976 _delta3>erase(e); 977 } 978 979 int vt = _tree_set>find(vb); 980 981 if (vt != tree) { 982 983 Edge r = _ugraph.oppositeEdge(e); 984 985 typename std::map<int, Edge>::iterator it = 986 (*_node_data)[ni].heap_index.find(vt); 987 988 if (it != (*_node_data)[ni].heap_index.end()) { 989 if ((*_node_data)[ni].heap[it>second] > rw) { 990 (*_node_data)[ni].heap.replace(it>second, r); 991 (*_node_data)[ni].heap.decrease(r, rw); 992 it>second = r; 993 } 994 } else { 995 (*_node_data)[ni].heap.push(r, rw); 996 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); 997 } 998 999 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { 1000 _blossom_set>decrease(n, (*_node_data)[ni].heap.prio()); 1001 1002 if (_delta2>state(blossom) != _delta2>IN_HEAP) { 1003 _delta2>push(blossom, _blossom_set>classPrio(blossom)  1004 (*_blossom_data)[blossom].offset); 1005 } else if ((*_delta2)[blossom] > 1006 _blossom_set>classPrio(blossom)  1007 (*_blossom_data)[blossom].offset){ 1008 _delta2>decrease(blossom, _blossom_set>classPrio(blossom)  1009 (*_blossom_data)[blossom].offset); 1010 } 1011 } 1012 } 1013 1014 } else if ((*_blossom_data)[vb].status == UNMATCHED) { 1015 if (_delta3>state(e) == _delta3>IN_HEAP) { 1016 _delta3>erase(e); 1017 } 1018 } else { 1019 1020 typename std::map<int, Edge>::iterator it = 1021 (*_node_data)[vi].heap_index.find(tree); 1022 1023 if (it != (*_node_data)[vi].heap_index.end()) { 1024 (*_node_data)[vi].heap.erase(it>second); 1025 (*_node_data)[vi].heap_index.erase(it); 1026 if ((*_node_data)[vi].heap.empty()) { 1027 _blossom_set>increase(v, std::numeric_limits<Value>::max()); 1028 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { 1029 _blossom_set>increase(v, (*_node_data)[vi].heap.prio()); 1030 } 1031 1032 if ((*_blossom_data)[vb].status == MATCHED) { 1033 if (_blossom_set>classPrio(vb) == 1034 std::numeric_limits<Value>::max()) { 1035 _delta2>erase(vb); 1036 } else if ((*_delta2)[vb] < _blossom_set>classPrio(vb)  1037 (*_blossom_data)[vb].offset) { 1038 _delta2>increase(vb, _blossom_set>classPrio(vb)  1039 (*_blossom_data)[vb].offset); 1040 } 1041 } 1042 } 1043 } 1044 } 1045 } 1046 } 1047 1048 void oddToMatched(int blossom) { 1049 (*_blossom_data)[blossom].offset = _delta_sum; 1050 1051 if (_blossom_set>classPrio(blossom) != 1052 std::numeric_limits<Value>::max()) { 1053 _delta2>push(blossom, _blossom_set>classPrio(blossom)  1054 (*_blossom_data)[blossom].offset); 1055 } 1056 1057 if (!_blossom_set>trivial(blossom)) { 1058 _delta4>erase(blossom); 1059 } 1060 } 1061 1062 void oddToEven(int blossom, int tree) { 1063 if (!_blossom_set>trivial(blossom)) { 1064 _delta4>erase(blossom); 1065 (*_blossom_data)[blossom].pot = 1066 2 * (2 * _delta_sum  (*_blossom_data)[blossom].offset); 1067 } 1068 1069 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 1070 n != INVALID; ++n) { 1071 int ni = (*_node_index)[n]; 1072 1073 _blossom_set>increase(n, std::numeric_limits<Value>::max()); 1074 1075 (*_node_data)[ni].heap.clear(); 1076 (*_node_data)[ni].heap_index.clear(); 1077 (*_node_data)[ni].pot += 1078 2 * _delta_sum  (*_blossom_data)[blossom].offset; 1079 1080 _delta1>push(n, (*_node_data)[ni].pot); 1081 1082 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { 1083 Node v = _ugraph.source(e); 1084 int vb = _blossom_set>find(v); 1085 int vi = (*_node_index)[v]; 1086 1087 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot  1088 dualScale * _weight[e]; 1089 1090 if ((*_blossom_data)[vb].status == EVEN) { 1091 if (_delta3>state(e) != _delta3>IN_HEAP && blossom != vb) { 1092 _delta3>push(e, rw / 2); 1093 } 1094 } else if ((*_blossom_data)[vb].status == UNMATCHED) { 1095 if (_delta3>state(e) != _delta3>IN_HEAP) { 1096 _delta3>push(e, rw); 1097 } 1098 } else { 1099 1100 typename std::map<int, Edge>::iterator it = 1101 (*_node_data)[vi].heap_index.find(tree); 1102 1103 if (it != (*_node_data)[vi].heap_index.end()) { 1104 if ((*_node_data)[vi].heap[it>second] > rw) { 1105 (*_node_data)[vi].heap.replace(it>second, e); 1106 (*_node_data)[vi].heap.decrease(e, rw); 1107 it>second = e; 1108 } 1109 } else { 1110 (*_node_data)[vi].heap.push(e, rw); 1111 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); 1112 } 1113 1114 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { 1115 _blossom_set>decrease(v, (*_node_data)[vi].heap.prio()); 1116 1117 if ((*_blossom_data)[vb].status == MATCHED) { 1118 if (_delta2>state(vb) != _delta2>IN_HEAP) { 1119 _delta2>push(vb, _blossom_set>classPrio(vb)  1120 (*_blossom_data)[vb].offset); 1121 } else if ((*_delta2)[vb] > _blossom_set>classPrio(vb)  1122 (*_blossom_data)[vb].offset) { 1123 _delta2>decrease(vb, _blossom_set>classPrio(vb)  1124 (*_blossom_data)[vb].offset); 1125 } 1126 } 1127 } 1128 } 1129 } 1130 } 1131 (*_blossom_data)[blossom].offset = 0; 1132 } 1133 1134 1135 void matchedToUnmatched(int blossom) { 1136 if (_delta2>state(blossom) == _delta2>IN_HEAP) { 1137 _delta2>erase(blossom); 1138 } 1139 1140 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 1141 n != INVALID; ++n) { 1142 int ni = (*_node_index)[n]; 1143 1144 _blossom_set>increase(n, std::numeric_limits<Value>::max()); 1145 1146 (*_node_data)[ni].heap.clear(); 1147 (*_node_data)[ni].heap_index.clear(); 1148 1149 for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) { 1150 Node v = _ugraph.target(e); 1151 int vb = _blossom_set>find(v); 1152 int vi = (*_node_index)[v]; 1153 1154 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot  1155 dualScale * _weight[e]; 1156 1157 if ((*_blossom_data)[vb].status == EVEN) { 1158 if (_delta3>state(e) != _delta3>IN_HEAP) { 1159 _delta3>push(e, rw); 1160 } 1161 } 1162 } 1163 } 1164 } 1165 1166 void unmatchedToMatched(int blossom) { 1167 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 1168 n != INVALID; ++n) { 1169 int ni = (*_node_index)[n]; 1170 1171 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { 1172 Node v = _ugraph.source(e); 1173 int vb = _blossom_set>find(v); 1174 int vi = (*_node_index)[v]; 1175 1176 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot  1177 dualScale * _weight[e]; 1178 1179 if (vb == blossom) { 1180 if (_delta3>state(e) == _delta3>IN_HEAP) { 1181 _delta3>erase(e); 1182 } 1183 } else if ((*_blossom_data)[vb].status == EVEN) { 1184 1185 if (_delta3>state(e) == _delta3>IN_HEAP) { 1186 _delta3>erase(e); 1187 } 1188 1189 int vt = _tree_set>find(vb); 1190 1191 Edge r = _ugraph.oppositeEdge(e); 1192 1193 typename std::map<int, Edge>::iterator it = 1194 (*_node_data)[ni].heap_index.find(vt); 1195 1196 if (it != (*_node_data)[ni].heap_index.end()) { 1197 if ((*_node_data)[ni].heap[it>second] > rw) { 1198 (*_node_data)[ni].heap.replace(it>second, r); 1199 (*_node_data)[ni].heap.decrease(r, rw); 1200 it>second = r; 1201 } 1202 } else { 1203 (*_node_data)[ni].heap.push(r, rw); 1204 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); 1205 } 1206 1207 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { 1208 _blossom_set>decrease(n, (*_node_data)[ni].heap.prio()); 1209 1210 if (_delta2>state(blossom) != _delta2>IN_HEAP) { 1211 _delta2>push(blossom, _blossom_set>classPrio(blossom)  1212 (*_blossom_data)[blossom].offset); 1213 } else if ((*_delta2)[blossom] > _blossom_set>classPrio(blossom) 1214 (*_blossom_data)[blossom].offset){ 1215 _delta2>decrease(blossom, _blossom_set>classPrio(blossom)  1216 (*_blossom_data)[blossom].offset); 1217 } 1218 } 1219 1220 } else if ((*_blossom_data)[vb].status == UNMATCHED) { 1221 if (_delta3>state(e) == _delta3>IN_HEAP) { 1222 _delta3>erase(e); 1223 } 1224 } 1225 } 1226 } 1227 } 1228 1229 void alternatePath(int even, int tree) { 1230 int odd; 1231 1232 evenToMatched(even, tree); 1233 (*_blossom_data)[even].status = MATCHED; 1234 1235 while ((*_blossom_data)[even].pred != INVALID) { 1236 odd = _blossom_set>find(_ugraph.target((*_blossom_data)[even].pred)); 1237 (*_blossom_data)[odd].status = MATCHED; 1238 oddToMatched(odd); 1239 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; 1240 1241 even = _blossom_set>find(_ugraph.target((*_blossom_data)[odd].pred)); 1242 (*_blossom_data)[even].status = MATCHED; 1243 evenToMatched(even, tree); 1244 (*_blossom_data)[even].next = 1245 _ugraph.oppositeEdge((*_blossom_data)[odd].pred); 1246 } 1247 1248 } 1249 1250 void destroyTree(int tree) { 1251 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { 1252 if ((*_blossom_data)[b].status == EVEN) { 1253 (*_blossom_data)[b].status = MATCHED; 1254 evenToMatched(b, tree); 1255 } else if ((*_blossom_data)[b].status == ODD) { 1256 (*_blossom_data)[b].status = MATCHED; 1257 oddToMatched(b); 1258 } 1259 } 1260 _tree_set>eraseClass(tree); 1261 } 1262 1263 1264 void unmatchNode(const Node& node) { 1265 int blossom = _blossom_set>find(node); 1266 int tree = _tree_set>find(blossom); 1267 1268 alternatePath(blossom, tree); 1269 destroyTree(tree); 1270 1271 (*_blossom_data)[blossom].status = UNMATCHED; 1272 (*_blossom_data)[blossom].base = node; 1273 matchedToUnmatched(blossom); 1274 } 1275 1276 1277 void augmentOnEdge(const UEdge& edge) { 1278 1279 int left = _blossom_set>find(_ugraph.source(edge)); 1280 int right = _blossom_set>find(_ugraph.target(edge)); 1281 1282 if ((*_blossom_data)[left].status == EVEN) { 1283 int left_tree = _tree_set>find(left); 1284 alternatePath(left, left_tree); 1285 destroyTree(left_tree); 1286 } else { 1287 (*_blossom_data)[left].status = MATCHED; 1288 unmatchedToMatched(left); 1289 } 1290 1291 if ((*_blossom_data)[right].status == EVEN) { 1292 int right_tree = _tree_set>find(right); 1293 alternatePath(right, right_tree); 1294 destroyTree(right_tree); 1295 } else { 1296 (*_blossom_data)[right].status = MATCHED; 1297 unmatchedToMatched(right); 1298 } 1299 1300 (*_blossom_data)[left].next = _ugraph.direct(edge, true); 1301 (*_blossom_data)[right].next = _ugraph.direct(edge, false); 1302 } 1303 1304 void extendOnEdge(const Edge& edge) { 1305 int base = _blossom_set>find(_ugraph.target(edge)); 1306 int tree = _tree_set>find(base); 1307 1308 int odd = _blossom_set>find(_ugraph.source(edge)); 1309 _tree_set>insert(odd, tree); 1310 (*_blossom_data)[odd].status = ODD; 1311 matchedToOdd(odd); 1312 (*_blossom_data)[odd].pred = edge; 1313 1314 int even = _blossom_set>find(_ugraph.target((*_blossom_data)[odd].next)); 1315 (*_blossom_data)[even].pred = (*_blossom_data)[even].next; 1316 _tree_set>insert(even, tree); 1317 (*_blossom_data)[even].status = EVEN; 1318 matchedToEven(even, tree); 1319 } 1320 1321 void shrinkOnEdge(const UEdge& uedge, int tree) { 1322 int nca = 1; 1323 std::vector<int> left_path, right_path; 1324 1325 { 1326 std::set<int> left_set, right_set; 1327 int left = _blossom_set>find(_ugraph.source(uedge)); 1328 left_path.push_back(left); 1329 left_set.insert(left); 1330 1331 int right = _blossom_set>find(_ugraph.target(uedge)); 1332 right_path.push_back(right); 1333 right_set.insert(right); 1334 1335 while (true) { 1336 1337 if ((*_blossom_data)[left].pred == INVALID) break; 1338 1339 left = 1340 _blossom_set>find(_ugraph.target((*_blossom_data)[left].pred)); 1341 left_path.push_back(left); 1342 left = 1343 _blossom_set>find(_ugraph.target((*_blossom_data)[left].pred)); 1344 left_path.push_back(left); 1345 1346 left_set.insert(left); 1347 1348 if (right_set.find(left) != right_set.end()) { 1349 nca = left; 1350 break; 1351 } 1352 1353 if ((*_blossom_data)[right].pred == INVALID) break; 1354 1355 right = 1356 _blossom_set>find(_ugraph.target((*_blossom_data)[right].pred)); 1357 right_path.push_back(right); 1358 right = 1359 _blossom_set>find(_ugraph.target((*_blossom_data)[right].pred)); 1360 right_path.push_back(right); 1361 1362 right_set.insert(right); 1363 1364 if (left_set.find(right) != left_set.end()) { 1365 nca = right; 1366 break; 1367 } 1368 1369 } 1370 1371 if (nca == 1) { 1372 if ((*_blossom_data)[left].pred == INVALID) { 1373 nca = right; 1374 while (left_set.find(nca) == left_set.end()) { 1375 nca = 1376 _blossom_set>find(_ugraph.target((*_blossom_data)[nca].pred)); 1377 right_path.push_back(nca); 1378 nca = 1379 _blossom_set>find(_ugraph.target((*_blossom_data)[nca].pred)); 1380 right_path.push_back(nca); 1381 } 1382 } else { 1383 nca = left; 1384 while (right_set.find(nca) == right_set.end()) { 1385 nca = 1386 _blossom_set>find(_ugraph.target((*_blossom_data)[nca].pred)); 1387 left_path.push_back(nca); 1388 nca = 1389 _blossom_set>find(_ugraph.target((*_blossom_data)[nca].pred)); 1390 left_path.push_back(nca); 1391 } 1392 } 1393 } 1394 } 1395 1396 std::vector<int> subblossoms; 1397 Edge prev; 1398 1399 prev = _ugraph.direct(uedge, true); 1400 for (int i = 0; left_path[i] != nca; i += 2) { 1401 subblossoms.push_back(left_path[i]); 1402 (*_blossom_data)[left_path[i]].next = prev; 1403 _tree_set>erase(left_path[i]); 1404 1405 subblossoms.push_back(left_path[i + 1]); 1406 (*_blossom_data)[left_path[i + 1]].status = EVEN; 1407 oddToEven(left_path[i + 1], tree); 1408 _tree_set>erase(left_path[i + 1]); 1409 prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred); 1410 } 1411 1412 int k = 0; 1413 while (right_path[k] != nca) ++k; 1414 1415 subblossoms.push_back(nca); 1416 (*_blossom_data)[nca].next = prev; 1417 1418 for (int i = k  2; i >= 0; i = 2) { 1419 subblossoms.push_back(right_path[i + 1]); 1420 (*_blossom_data)[right_path[i + 1]].status = EVEN; 1421 oddToEven(right_path[i + 1], tree); 1422 _tree_set>erase(right_path[i + 1]); 1423 1424 (*_blossom_data)[right_path[i + 1]].next = 1425 (*_blossom_data)[right_path[i + 1]].pred; 1426 1427 subblossoms.push_back(right_path[i]); 1428 _tree_set>erase(right_path[i]); 1429 } 1430 1431 int surface = 1432 _blossom_set>join(subblossoms.begin(), subblossoms.end()); 1433 1434 for (int i = 0; i < int(subblossoms.size()); ++i) { 1435 if (!_blossom_set>trivial(subblossoms[i])) { 1436 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; 1437 } 1438 (*_blossom_data)[subblossoms[i]].status = MATCHED; 1439 } 1440 1441 (*_blossom_data)[surface].pot = 2 * _delta_sum; 1442 (*_blossom_data)[surface].offset = 0; 1443 (*_blossom_data)[surface].status = EVEN; 1444 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; 1445 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; 1446 1447 _tree_set>insert(surface, tree); 1448 _tree_set>erase(nca); 1449 } 1450 1451 void splitBlossom(int blossom) { 1452 Edge next = (*_blossom_data)[blossom].next; 1453 Edge pred = (*_blossom_data)[blossom].pred; 1454 1455 int tree = _tree_set>find(blossom); 1456 1457 (*_blossom_data)[blossom].status = MATCHED; 1458 oddToMatched(blossom); 1459 if (_delta2>state(blossom) == _delta2>IN_HEAP) { 1460 _delta2>erase(blossom); 1461 } 1462 1463 std::vector<int> subblossoms; 1464 _blossom_set>split(blossom, std::back_inserter(subblossoms)); 1465 1466 Value offset = (*_blossom_data)[blossom].offset; 1467 int b = _blossom_set>find(_ugraph.source(pred)); 1468 int d = _blossom_set>find(_ugraph.source(next)); 1469 1470 int ib, id; 1471 for (int i = 0; i < int(subblossoms.size()); ++i) { 1472 if (subblossoms[i] == b) ib = i; 1473 if (subblossoms[i] == d) id = i; 1474 1475 (*_blossom_data)[subblossoms[i]].offset = offset; 1476 if (!_blossom_set>trivial(subblossoms[i])) { 1477 (*_blossom_data)[subblossoms[i]].pot = 2 * offset; 1478 } 1479 if (_blossom_set>classPrio(subblossoms[i]) != 1480 std::numeric_limits<Value>::max()) { 1481 _delta2>push(subblossoms[i], 1482 _blossom_set>classPrio(subblossoms[i])  1483 (*_blossom_data)[subblossoms[i]].offset); 1484 } 1485 } 1486 1487 if (id > ib ? ((id  ib) % 2 == 0) : ((ib  id) % 2 == 1)) { 1488 for (int i = (id + 1) % subblossoms.size(); 1489 i != ib; i = (i + 2) % subblossoms.size()) { 1490 int sb = subblossoms[i]; 1491 int tb = subblossoms[(i + 1) % subblossoms.size()]; 1492 (*_blossom_data)[sb].next = 1493 _ugraph.oppositeEdge((*_blossom_data)[tb].next); 1494 } 1495 1496 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { 1497 int sb = subblossoms[i]; 1498 int tb = subblossoms[(i + 1) % subblossoms.size()]; 1499 int ub = subblossoms[(i + 2) % subblossoms.size()]; 1500 1501 (*_blossom_data)[sb].status = ODD; 1502 matchedToOdd(sb); 1503 _tree_set>insert(sb, tree); 1504 (*_blossom_data)[sb].pred = pred; 1505 (*_blossom_data)[sb].next = 1506 _ugraph.oppositeEdge((*_blossom_data)[tb].next); 1507 1508 pred = (*_blossom_data)[ub].next; 1509 1510 (*_blossom_data)[tb].status = EVEN; 1511 matchedToEven(tb, tree); 1512 _tree_set>insert(tb, tree); 1513 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; 1514 } 1515 1516 (*_blossom_data)[subblossoms[id]].status = ODD; 1517 matchedToOdd(subblossoms[id]); 1518 _tree_set>insert(subblossoms[id], tree); 1519 (*_blossom_data)[subblossoms[id]].next = next; 1520 (*_blossom_data)[subblossoms[id]].pred = pred; 1521 1522 } else { 1523 1524 for (int i = (ib + 1) % subblossoms.size(); 1525 i != id; i = (i + 2) % subblossoms.size()) { 1526 int sb = subblossoms[i]; 1527 int tb = subblossoms[(i + 1) % subblossoms.size()]; 1528 (*_blossom_data)[sb].next = 1529 _ugraph.oppositeEdge((*_blossom_data)[tb].next); 1530 } 1531 1532 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { 1533 int sb = subblossoms[i]; 1534 int tb = subblossoms[(i + 1) % subblossoms.size()]; 1535 int ub = subblossoms[(i + 2) % subblossoms.size()]; 1536 1537 (*_blossom_data)[sb].status = ODD; 1538 matchedToOdd(sb); 1539 _tree_set>insert(sb, tree); 1540 (*_blossom_data)[sb].next = next; 1541 (*_blossom_data)[sb].pred = 1542 _ugraph.oppositeEdge((*_blossom_data)[tb].next); 1543 1544 (*_blossom_data)[tb].status = EVEN; 1545 matchedToEven(tb, tree); 1546 _tree_set>insert(tb, tree); 1547 (*_blossom_data)[tb].pred = 1548 (*_blossom_data)[tb].next = 1549 _ugraph.oppositeEdge((*_blossom_data)[ub].next); 1550 next = (*_blossom_data)[ub].next; 1551 } 1552 1553 (*_blossom_data)[subblossoms[ib]].status = ODD; 1554 matchedToOdd(subblossoms[ib]); 1555 _tree_set>insert(subblossoms[ib], tree); 1556 (*_blossom_data)[subblossoms[ib]].next = next; 1557 (*_blossom_data)[subblossoms[ib]].pred = pred; 1558 } 1559 _tree_set>erase(blossom); 1560 } 1561 1562 void extractBlossom(int blossom, const Node& base, const Edge& matching) { 1563 if (_blossom_set>trivial(blossom)) { 1564 int bi = (*_node_index)[base]; 1565 Value pot = (*_node_data)[bi].pot; 1566 1567 _matching>set(base, matching); 1568 _blossom_node_list.push_back(base); 1569 _node_potential>set(base, pot); 1570 } else { 1571 1572 Value pot = (*_blossom_data)[blossom].pot; 1573 int bn = _blossom_node_list.size(); 1574 1575 std::vector<int> subblossoms; 1576 _blossom_set>split(blossom, std::back_inserter(subblossoms)); 1577 int b = _blossom_set>find(base); 1578 int ib = 1; 1579 for (int i = 0; i < int(subblossoms.size()); ++i) { 1580 if (subblossoms[i] == b) { ib = i; break; } 1581 } 1582 1583 for (int i = 1; i < int(subblossoms.size()); i += 2) { 1584 int sb = subblossoms[(ib + i) % subblossoms.size()]; 1585 int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; 1586 1587 Edge m = (*_blossom_data)[tb].next; 1588 extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m)); 1589 extractBlossom(tb, _ugraph.source(m), m); 1590 } 1591 extractBlossom(subblossoms[ib], base, matching); 1592 1593 int en = _blossom_node_list.size(); 1594 1595 _blossom_potential.push_back(BlossomVariable(bn, en, pot)); 1596 } 1597 } 1598 1599 void extractMatching() { 1600 std::vector<int> blossoms; 1601 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { 1602 blossoms.push_back(c); 1603 } 1604 1605 for (int i = 0; i < int(blossoms.size()); ++i) { 1606 if ((*_blossom_data)[blossoms[i]].status == MATCHED) { 1607 1608 Value offset = (*_blossom_data)[blossoms[i]].offset; 1609 (*_blossom_data)[blossoms[i]].pot += 2 * offset; 1610 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); 1611 n != INVALID; ++n) { 1612 (*_node_data)[(*_node_index)[n]].pot = offset; 1613 } 1614 1615 Edge matching = (*_blossom_data)[blossoms[i]].next; 1616 Node base = _ugraph.source(matching); 1617 extractBlossom(blossoms[i], base, matching); 1618 } else { 1619 Node base = (*_blossom_data)[blossoms[i]].base; 1620 extractBlossom(blossoms[i], base, INVALID); 1621 } 1622 } 1623 } 1624 1625 public: 1626 1627 /// \brief Constructor 1628 /// 1629 /// Constructor. 1630 MaxWeightedMatching(const UGraph& ugraph, const WeightMap& weight) 1631 : _ugraph(ugraph), _weight(weight), _matching(0), 1632 _node_potential(0), _blossom_potential(), _blossom_node_list(), 1633 _node_num(0), _blossom_num(0), 1634 1635 _blossom_index(0), _blossom_set(0), _blossom_data(0), 1636 _node_index(0), _node_heap_index(0), _node_data(0), 1637 _tree_set_index(0), _tree_set(0), 1638 1639 _delta1_index(0), _delta1(0), 1640 _delta2_index(0), _delta2(0), 1641 _delta3_index(0), _delta3(0), 1642 _delta4_index(0), _delta4(0), 1643 1644 _delta_sum() {} 1645 1646 ~MaxWeightedMatching() { 1647 destroyStructures(); 1648 } 1649 1650 /// \name Execution control 1651 /// The simplest way to execute the algorithm is to use the member 1652 /// \c run() member function. 1653 1654 ///@{ 1655 1656 /// \brief Initialize the algorithm 1657 /// 1658 /// Initialize the algorithm 1659 void init() { 1660 createStructures(); 1661 1662 for (EdgeIt e(_ugraph); e != INVALID; ++e) { 1663 _node_heap_index>set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP); 1664 } 1665 for (NodeIt n(_ugraph); n != INVALID; ++n) { 1666 _delta1_index>set(n, _delta1>PRE_HEAP); 1667 } 1668 for (UEdgeIt e(_ugraph); e != INVALID; ++e) { 1669 _delta3_index>set(e, _delta3>PRE_HEAP); 1670 } 1671 for (int i = 0; i < _blossom_num; ++i) { 1672 _delta2_index>set(i, _delta2>PRE_HEAP); 1673 _delta4_index>set(i, _delta4>PRE_HEAP); 1674 } 1675 1676 int index = 0; 1677 for (NodeIt n(_ugraph); n != INVALID; ++n) { 1678 Value max = 0; 1679 for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) { 1680 if (_ugraph.target(e) == n) continue; 1681 if ((dualScale * _weight[e]) / 2 > max) { 1682 max = (dualScale * _weight[e]) / 2; 1683 } 1684 } 1685 _node_index>set(n, index); 1686 (*_node_data)[index].pot = max; 1687 _delta1>push(n, max); 1688 int blossom = 1689 _blossom_set>insert(n, std::numeric_limits<Value>::max()); 1690 1691 _tree_set>insert(blossom); 1692 1693 (*_blossom_data)[blossom].status = EVEN; 1694 (*_blossom_data)[blossom].pred = INVALID; 1695 (*_blossom_data)[blossom].next = INVALID; 1696 (*_blossom_data)[blossom].pot = 0; 1697 (*_blossom_data)[blossom].offset = 0; 1698 ++index; 1699 } 1700 for (UEdgeIt e(_ugraph); e != INVALID; ++e) { 1701 int si = (*_node_index)[_ugraph.source(e)]; 1702 int ti = (*_node_index)[_ugraph.target(e)]; 1703 if (_ugraph.source(e) != _ugraph.target(e)) { 1704 _delta3>push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot  1705 dualScale * _weight[e]) / 2); 1706 } 1707 } 1708 } 1709 1710 /// \brief Starts the algorithm 1711 /// 1712 /// Starts the algorithm 1713 void start() { 1714 enum OpType { 1715 D1, D2, D3, D4 1716 }; 1717 1718 int unmatched = _node_num; 1719 while (unmatched > 0) { 1720 Value d1 = !_delta1>empty() ? 1721 _delta1>prio() : std::numeric_limits<Value>::max(); 1722 1723 Value d2 = !_delta2>empty() ? 1724 _delta2>prio() : std::numeric_limits<Value>::max(); 1725 1726 Value d3 = !_delta3>empty() ? 1727 _delta3>prio() : std::numeric_limits<Value>::max(); 1728 1729 Value d4 = !_delta4>empty() ? 1730 _delta4>prio() : std::numeric_limits<Value>::max(); 1731 1732 _delta_sum = d1; OpType ot = D1; 1733 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } 1734 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } 1735 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } 1736 1737 1738 switch (ot) { 1739 case D1: 1740 { 1741 Node n = _delta1>top(); 1742 unmatchNode(n); 1743 unmatched; 1744 } 1745 break; 1746 case D2: 1747 { 1748 int blossom = _delta2>top(); 1749 Node n = _blossom_set>classTop(blossom); 1750 Edge e = (*_node_data)[(*_node_index)[n]].heap.top(); 1751 extendOnEdge(e); 1752 } 1753 break; 1754 case D3: 1755 { 1756 UEdge e = _delta3>top(); 1757 1758 int left_blossom = _blossom_set>find(_ugraph.source(e)); 1759 int right_blossom = _blossom_set>find(_ugraph.target(e)); 1760 1761 if (left_blossom == right_blossom) { 1762 _delta3>pop(); 1763 } else { 1764 int left_tree; 1765 if ((*_blossom_data)[left_blossom].status == EVEN) { 1766 left_tree = _tree_set>find(left_blossom); 1767 } else { 1768 left_tree = 1; 1769 ++unmatched; 1770 } 1771 int right_tree; 1772 if ((*_blossom_data)[right_blossom].status == EVEN) { 1773 right_tree = _tree_set>find(right_blossom); 1774 } else { 1775 right_tree = 1; 1776 ++unmatched; 1777 } 1778 1779 if (left_tree == right_tree) { 1780 shrinkOnEdge(e, left_tree); 1781 } else { 1782 augmentOnEdge(e); 1783 unmatched = 2; 1784 } 1785 } 1786 } break; 1787 case D4: 1788 splitBlossom(_delta4>top()); 1789 break; 1790 } 1791 } 1792 extractMatching(); 1793 } 1794 1795 /// \brief Runs %MaxWeightedMatching algorithm. 1796 /// 1797 /// This method runs the %MaxWeightedMatching algorithm. 1798 /// 1799 /// \note mwm.run() is just a shortcut of the following code. 1800 /// \code 1801 /// mwm.init(); 1802 /// mwm.start(); 1803 /// \endcode 1804 void run() { 1805 init(); 1806 start(); 1807 } 1808 1809 /// @} 1810 1811 /// \name Primal solution 1812 /// Functions for get the primal solution, ie. the matching. 1813 1814 /// @{ 1815 1816 /// \brief Returns the matching value. 1817 /// 1818 /// Returns the matching value. 1819 Value matchingValue() const { 1820 Value sum = 0; 1821 for (NodeIt n(_ugraph); n != INVALID; ++n) { 1822 if ((*_matching)[n] != INVALID) { 1823 sum += _weight[(*_matching)[n]]; 1824 } 1825 } 1826 return sum /= 2; 1827 } 1828 1829 /// \brief Returns true when the edge is in the matching. 1830 /// 1831 /// Returns true when the edge is in the matching. 1832 bool matching(const UEdge& edge) const { 1833 return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true); 1834 } 1835 1836 /// \brief Returns the incident matching edge. 1837 /// 1838 /// Returns the incident matching edge from given node. If the 1839 /// node is not matched then it gives back \c INVALID. 1840 Edge matching(const Node& node) const { 1841 return (*_matching)[node]; 1842 } 1843 1844 /// \brief Returns the mate of the node. 1845 /// 1846 /// Returns the adjancent node in a mathcing edge. If the node is 1847 /// not matched then it gives back \c INVALID. 1848 Node mate(const Node& node) const { 1849 return (*_matching)[node] != INVALID ? 1850 _ugraph.target((*_matching)[node]) : INVALID; 1851 } 1852 1853 /// @} 1854 1855 /// \name Dual solution 1856 /// Functions for get the dual solution. 1857 1858 /// @{ 1859 1860 /// \brief Returns the value of the dual solution. 1861 /// 1862 /// Returns the value of the dual solution. It should be equal to 1863 /// the primal value scaled by \ref dualScale "dual scale". 1864 Value dualValue() const { 1865 Value sum = 0; 1866 for (NodeIt n(_ugraph); n != INVALID; ++n) { 1867 sum += nodeValue(n); 1868 } 1869 for (int i = 0; i < blossomNum(); ++i) { 1870 sum += blossomValue(i) * (blossomSize(i) / 2); 1871 } 1872 return sum; 1873 } 1874 1875 /// \brief Returns the value of the node. 1876 /// 1877 /// Returns the the value of the node. 1878 Value nodeValue(const Node& n) const { 1879 return (*_node_potential)[n]; 1880 } 1881 1882 /// \brief Returns the number of the blossoms in the basis. 1883 /// 1884 /// Returns the number of the blossoms in the basis. 1885 /// \see BlossomIt 1886 int blossomNum() const { 1887 return _blossom_potential.size(); 1888 } 1889 1890 1891 /// \brief Returns the number of the nodes in the blossom. 1892 /// 1893 /// Returns the number of the nodes in the blossom. 1894 int blossomSize(int k) const { 1895 return _blossom_potential[k].end  _blossom_potential[k].begin; 1896 } 1897 1898 /// \brief Returns the value of the blossom. 1899 /// 1900 /// Returns the the value of the blossom. 1901 /// \see BlossomIt 1902 Value blossomValue(int k) const { 1903 return _blossom_potential[k].value; 1904 } 1905 1906 /// \brief Lemon iterator for get the items of the blossom. 1907 /// 1908 /// Lemon iterator for get the nodes of the blossom. This class 1909 /// provides a common style lemon iterator which gives back a 1910 /// subset of the nodes. 1911 class BlossomIt { 1912 public: 1913 1914 /// \brief Constructor. 1915 /// 1916 /// Constructor for get the nodes of the variable. 1917 BlossomIt(const MaxWeightedMatching& algorithm, int variable) 1918 : _algorithm(&algorithm) 1919 { 1920 _index = _algorithm>_blossom_potential[variable].begin; 1921 _last = _algorithm>_blossom_potential[variable].end; 1922 } 1923 1924 /// \brief Invalid constructor. 1925 /// 1926 /// Invalid constructor. 1927 BlossomIt(Invalid) : _index(1) {} 1928 1929 /// \brief Conversion to node. 1930 /// 1931 /// Conversion to node. 1932 operator Node() const { 1933 return _algorithm ? _algorithm>_blossom_node_list[_index] : INVALID; 1934 } 1935 1936 /// \brief Increment operator. 1937 /// 1938 /// Increment operator. 1939 BlossomIt& operator++() { 1940 ++_index; 1941 if (_index == _last) { 1942 _index = 1; 1943 } 1944 return *this; 1945 } 1946 1947 bool operator==(const BlossomIt& it) const { 1948 return _index == it._index; 1949 } 1950 bool operator!=(const BlossomIt& it) const { 1951 return _index != it._index; 1952 } 1953 1954 private: 1955 const MaxWeightedMatching* _algorithm; 1956 int _last; 1957 int _index; 1958 }; 1959 1960 /// @} 1961 1962 }; 1963 1964 /// \ingroup matching 1965 /// 1966 /// \brief Weighted perfect matching in general undirected graphs 1967 /// 1968 /// This class provides an efficient implementation of Edmond's 1969 /// maximum weighted perfecr matching algorithm. The implementation 1970 /// is based on extensive use of priority queues and provides 1971 /// \f$O(nm\log(n))\f$ time complexity. 1972 /// 1973 /// The maximum weighted matching problem is to find undirected 1974 /// edges in the graph with maximum overall weight and no two of 1975 /// them shares their endpoints and covers all nodes. The problem 1976 /// can be formulated with the next linear program: 1977 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f] 1978 ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert  1}{2} \quad \forall B\in\mathcal{O}\f] 1979 /// \f[x_e \ge 0\quad \forall e\in E\f] 1980 /// \f[\max \sum_{e\in E}x_ew_e\f] 1981 /// where \f$\delta(X)\f$ is the set of edges incident to a node in 1982 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in 1983 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of 1984 /// the nodes. 1985 /// 1986 /// The algorithm calculates an optimal matching and a proof of the 1987 /// optimality. The solution of the dual problem can be used to check 1988 /// the result of the algorithm. The dual linear problem is the next: 1989 /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f] 1990 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] 1991 /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert  1}{2}z_B\f] 1992 /// 1993 /// The algorithm can be executed with \c run() or the \c init() and 1994 /// then the \c start() member functions. After it the matching can 1995 /// be asked with \c matching() or mate() functions. The dual 1996 /// solution can be get with \c nodeValue(), \c blossomNum() and \c 1997 /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt 1998 /// "BlossomIt" nested class which is able to iterate on the nodes 1999 /// of a blossom. If the value type is integral then the dual 2000 /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4". 2001 /// 2002 /// \author Balazs Dezso 2003 template <typename _UGraph, 2004 typename _WeightMap = typename _UGraph::template UEdgeMap<int> > 2005 class MaxWeightedPerfectMatching { 2006 public: 2007 2008 typedef _UGraph UGraph; 2009 typedef _WeightMap WeightMap; 2010 typedef typename WeightMap::Value Value; 2011 2012 /// \brief Scaling factor for dual solution 2013 /// 2014 /// Scaling factor for dual solution, it is equal to 4 or 1 2015 /// according to the value type. 2016 static const int dualScale = 2017 std::numeric_limits<Value>::is_integer ? 4 : 1; 2018 2019 typedef typename UGraph::template NodeMap<typename UGraph::Edge> 2020 MatchingMap; 2021 2022 private: 2023 2024 UGRAPH_TYPEDEFS(typename UGraph); 2025 2026 typedef typename UGraph::template NodeMap<Value> NodePotential; 2027 typedef std::vector<Node> BlossomNodeList; 2028 2029 struct BlossomVariable { 2030 int begin, end; 2031 Value value; 2032 2033 BlossomVariable(int _begin, int _end, Value _value) 2034 : begin(_begin), end(_end), value(_value) {} 2035 2036 }; 2037 2038 typedef std::vector<BlossomVariable> BlossomPotential; 2039 2040 const UGraph& _ugraph; 2041 const WeightMap& _weight; 2042 2043 MatchingMap* _matching; 2044 2045 NodePotential* _node_potential; 2046 2047 BlossomPotential _blossom_potential; 2048 BlossomNodeList _blossom_node_list; 2049 2050 int _node_num; 2051 int _blossom_num; 2052 2053 typedef typename UGraph::template NodeMap<int> NodeIntMap; 2054 typedef typename UGraph::template EdgeMap<int> EdgeIntMap; 2055 typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap; 2056 typedef IntegerMap<int> IntIntMap; 2057 2058 enum Status { 2059 EVEN = 1, MATCHED = 0, ODD = 1 2060 }; 2061 2062 typedef HeapUnionFind<Value, NodeIntMap> BlossomSet; 2063 struct BlossomData { 2064 int tree; 2065 Status status; 2066 Edge pred, next; 2067 Value pot, offset; 2068 }; 2069 2070 NodeIntMap *_blossom_index; 2071 BlossomSet *_blossom_set; 2072 IntegerMap<BlossomData>* _blossom_data; 2073 2074 NodeIntMap *_node_index; 2075 EdgeIntMap *_node_heap_index; 2076 2077 struct NodeData { 2078 2079 NodeData(EdgeIntMap& node_heap_index) 2080 : heap(node_heap_index) {} 2081 2082 int blossom; 2083 Value pot; 2084 BinHeap<Value, EdgeIntMap> heap; 2085 std::map<int, Edge> heap_index; 2086 2087 int tree; 2088 }; 2089 2090 IntegerMap<NodeData>* _node_data; 2091 2092 typedef ExtendFindEnum<IntIntMap> TreeSet; 2093 2094 IntIntMap *_tree_set_index; 2095 TreeSet *_tree_set; 2096 2097 IntIntMap *_delta2_index; 2098 BinHeap<Value, IntIntMap> *_delta2; 2099 2100 UEdgeIntMap *_delta3_index; 2101 BinHeap<Value, UEdgeIntMap> *_delta3; 2102 2103 IntIntMap *_delta4_index; 2104 BinHeap<Value, IntIntMap> *_delta4; 2105 2106 Value _delta_sum; 2107 2108 void createStructures() { 2109 _node_num = countNodes(_ugraph); 2110 _blossom_num = _node_num * 3 / 2; 2111 2112 if (!_matching) { 2113 _matching = new MatchingMap(_ugraph); 2114 } 2115 if (!_node_potential) { 2116 _node_potential = new NodePotential(_ugraph); 2117 } 2118 if (!_blossom_set) { 2119 _blossom_index = new NodeIntMap(_ugraph); 2120 _blossom_set = new BlossomSet(*_blossom_index); 2121 _blossom_data = new IntegerMap<BlossomData>(_blossom_num); 2122 } 2123 2124 if (!_node_index) { 2125 _node_index = new NodeIntMap(_ugraph); 2126 _node_heap_index = new EdgeIntMap(_ugraph); 2127 _node_data = new IntegerMap<NodeData>(_node_num, 2128 NodeData(*_node_heap_index)); 2129 } 2130 2131 if (!_tree_set) { 2132 _tree_set_index = new IntIntMap(_blossom_num); 2133 _tree_set = new TreeSet(*_tree_set_index); 2134 } 2135 if (!_delta2) { 2136 _delta2_index = new IntIntMap(_blossom_num); 2137 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index); 2138 } 2139 if (!_delta3) { 2140 _delta3_index = new UEdgeIntMap(_ugraph); 2141 _delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index); 2142 } 2143 if (!_delta4) { 2144 _delta4_index = new IntIntMap(_blossom_num); 2145 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index); 2146 } 2147 } 2148 2149 void destroyStructures() { 2150 _node_num = countNodes(_ugraph); 2151 _blossom_num = _node_num * 3 / 2; 2152 2153 if (_matching) { 2154 delete _matching; 2155 } 2156 if (_node_potential) { 2157 delete _node_potential; 2158 } 2159 if (_blossom_set) { 2160 delete _blossom_index; 2161 delete _blossom_set; 2162 delete _blossom_data; 2163 } 2164 2165 if (_node_index) { 2166 delete _node_index; 2167 delete _node_heap_index; 2168 delete _node_data; 2169 } 2170 2171 if (_tree_set) { 2172 delete _tree_set_index; 2173 delete _tree_set; 2174 } 2175 if (_delta2) { 2176 delete _delta2_index; 2177 delete _delta2; 2178 } 2179 if (_delta3) { 2180 delete _delta3_index; 2181 delete _delta3; 2182 } 2183 if (_delta4) { 2184 delete _delta4_index; 2185 delete _delta4; 2186 } 2187 } 2188 2189 void matchedToEven(int blossom, int tree) { 2190 if (_delta2>state(blossom) == _delta2>IN_HEAP) { 2191 _delta2>erase(blossom); 2192 } 2193 2194 if (!_blossom_set>trivial(blossom)) { 2195 (*_blossom_data)[blossom].pot = 2196 2 * (_delta_sum  (*_blossom_data)[blossom].offset); 2197 } 2198 2199 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 2200 n != INVALID; ++n) { 2201 2202 _blossom_set>increase(n, std::numeric_limits<Value>::max()); 2203 int ni = (*_node_index)[n]; 2204 2205 (*_node_data)[ni].heap.clear(); 2206 (*_node_data)[ni].heap_index.clear(); 2207 2208 (*_node_data)[ni].pot += _delta_sum  (*_blossom_data)[blossom].offset; 2209 2210 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { 2211 Node v = _ugraph.source(e); 2212 int vb = _blossom_set>find(v); 2213 int vi = (*_node_index)[v]; 2214 2215 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot  2216 dualScale * _weight[e]; 2217 2218 if ((*_blossom_data)[vb].status == EVEN) { 2219 if (_delta3>state(e) != _delta3>IN_HEAP && blossom != vb) { 2220 _delta3>push(e, rw / 2); 2221 } 2222 } else { 2223 typename std::map<int, Edge>::iterator it = 2224 (*_node_data)[vi].heap_index.find(tree); 2225 2226 if (it != (*_node_data)[vi].heap_index.end()) { 2227 if ((*_node_data)[vi].heap[it>second] > rw) { 2228 (*_node_data)[vi].heap.replace(it>second, e); 2229 (*_node_data)[vi].heap.decrease(e, rw); 2230 it>second = e; 2231 } 2232 } else { 2233 (*_node_data)[vi].heap.push(e, rw); 2234 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); 2235 } 2236 2237 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { 2238 _blossom_set>decrease(v, (*_node_data)[vi].heap.prio()); 2239 2240 if ((*_blossom_data)[vb].status == MATCHED) { 2241 if (_delta2>state(vb) != _delta2>IN_HEAP) { 2242 _delta2>push(vb, _blossom_set>classPrio(vb)  2243 (*_blossom_data)[vb].offset); 2244 } else if ((*_delta2)[vb] > _blossom_set>classPrio(vb)  2245 (*_blossom_data)[vb].offset){ 2246 _delta2>decrease(vb, _blossom_set>classPrio(vb)  2247 (*_blossom_data)[vb].offset); 2248 } 2249 } 2250 } 2251 } 2252 } 2253 } 2254 (*_blossom_data)[blossom].offset = 0; 2255 } 2256 2257 void matchedToOdd(int blossom) { 2258 if (_delta2>state(blossom) == _delta2>IN_HEAP) { 2259 _delta2>erase(blossom); 2260 } 2261 (*_blossom_data)[blossom].offset += _delta_sum; 2262 if (!_blossom_set>trivial(blossom)) { 2263 _delta4>push(blossom, (*_blossom_data)[blossom].pot / 2 + 2264 (*_blossom_data)[blossom].offset); 2265 } 2266 } 2267 2268 void evenToMatched(int blossom, int tree) { 2269 if (!_blossom_set>trivial(blossom)) { 2270 (*_blossom_data)[blossom].pot += 2 * _delta_sum; 2271 } 2272 2273 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 2274 n != INVALID; ++n) { 2275 int ni = (*_node_index)[n]; 2276 (*_node_data)[ni].pot = _delta_sum; 2277 2278 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { 2279 Node v = _ugraph.source(e); 2280 int vb = _blossom_set>find(v); 2281 int vi = (*_node_index)[v]; 2282 2283 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot  2284 dualScale * _weight[e]; 2285 2286 if (vb == blossom) { 2287 if (_delta3>state(e) == _delta3>IN_HEAP) { 2288 _delta3>erase(e); 2289 } 2290 } else if ((*_blossom_data)[vb].status == EVEN) { 2291 2292 if (_delta3>state(e) == _delta3>IN_HEAP) { 2293 _delta3>erase(e); 2294 } 2295 2296 int vt = _tree_set>find(vb); 2297 2298 if (vt != tree) { 2299 2300 Edge r = _ugraph.oppositeEdge(e); 2301 2302 typename std::map<int, Edge>::iterator it = 2303 (*_node_data)[ni].heap_index.find(vt); 2304 2305 if (it != (*_node_data)[ni].heap_index.end()) { 2306 if ((*_node_data)[ni].heap[it>second] > rw) { 2307 (*_node_data)[ni].heap.replace(it>second, r); 2308 (*_node_data)[ni].heap.decrease(r, rw); 2309 it>second = r; 2310 } 2311 } else { 2312 (*_node_data)[ni].heap.push(r, rw); 2313 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); 2314 } 2315 2316 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { 2317 _blossom_set>decrease(n, (*_node_data)[ni].heap.prio()); 2318 2319 if (_delta2>state(blossom) != _delta2>IN_HEAP) { 2320 _delta2>push(blossom, _blossom_set>classPrio(blossom)  2321 (*_blossom_data)[blossom].offset); 2322 } else if ((*_delta2)[blossom] > 2323 _blossom_set>classPrio(blossom)  2324 (*_blossom_data)[blossom].offset){ 2325 _delta2>decrease(blossom, _blossom_set>classPrio(blossom)  2326 (*_blossom_data)[blossom].offset); 2327 } 2328 } 2329 } 2330 } else { 2331 2332 typename std::map<int, Edge>::iterator it = 2333 (*_node_data)[vi].heap_index.find(tree); 2334 2335 if (it != (*_node_data)[vi].heap_index.end()) { 2336 (*_node_data)[vi].heap.erase(it>second); 2337 (*_node_data)[vi].heap_index.erase(it); 2338 if ((*_node_data)[vi].heap.empty()) { 2339 _blossom_set>increase(v, std::numeric_limits<Value>::max()); 2340 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { 2341 _blossom_set>increase(v, (*_node_data)[vi].heap.prio()); 2342 } 2343 2344 if ((*_blossom_data)[vb].status == MATCHED) { 2345 if (_blossom_set>classPrio(vb) == 2346 std::numeric_limits<Value>::max()) { 2347 _delta2>erase(vb); 2348 } else if ((*_delta2)[vb] < _blossom_set>classPrio(vb)  2349 (*_blossom_data)[vb].offset) { 2350 _delta2>increase(vb, _blossom_set>classPrio(vb)  2351 (*_blossom_data)[vb].offset); 2352 } 2353 } 2354 } 2355 } 2356 } 2357 } 2358 } 2359 2360 void oddToMatched(int blossom) { 2361 (*_blossom_data)[blossom].offset = _delta_sum; 2362 2363 if (_blossom_set>classPrio(blossom) != 2364 std::numeric_limits<Value>::max()) { 2365 _delta2>push(blossom, _blossom_set>classPrio(blossom)  2366 (*_blossom_data)[blossom].offset); 2367 } 2368 2369 if (!_blossom_set>trivial(blossom)) { 2370 _delta4>erase(blossom); 2371 } 2372 } 2373 2374 void oddToEven(int blossom, int tree) { 2375 if (!_blossom_set>trivial(blossom)) { 2376 _delta4>erase(blossom); 2377 (*_blossom_data)[blossom].pot = 2378 2 * (2 * _delta_sum  (*_blossom_data)[blossom].offset); 2379 } 2380 2381 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 2382 n != INVALID; ++n) { 2383 int ni = (*_node_index)[n]; 2384 2385 _blossom_set>increase(n, std::numeric_limits<Value>::max()); 2386 2387 (*_node_data)[ni].heap.clear(); 2388 (*_node_data)[ni].heap_index.clear(); 2389 (*_node_data)[ni].pot += 2390 2 * _delta_sum  (*_blossom_data)[blossom].offset; 2391 2392 for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { 2393 Node v = _ugraph.source(e); 2394 int vb = _blossom_set>find(v); 2395 int vi = (*_node_index)[v]; 2396 2397 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot  2398 dualScale * _weight[e]; 2399 2400 if ((*_blossom_data)[vb].status == EVEN) { 2401 if (_delta3>state(e) != _delta3>IN_HEAP && blossom != vb) { 2402 _delta3>push(e, rw / 2); 2403 } 2404 } else { 2405 2406 typename std::map<int, Edge>::iterator it = 2407 (*_node_data)[vi].heap_index.find(tree); 2408 2409 if (it != (*_node_data)[vi].heap_index.end()) { 2410 if ((*_node_data)[vi].heap[it>second] > rw) { 2411 (*_node_data)[vi].heap.replace(it>second, e); 2412 (*_node_data)[vi].heap.decrease(e, rw); 2413 it>second = e; 2414 } 2415 } else { 2416 (*_node_data)[vi].heap.push(e, rw); 2417 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); 2418 } 2419 2420 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { 2421 _blossom_set>decrease(v, (*_node_data)[vi].heap.prio()); 2422 2423 if ((*_blossom_data)[vb].status == MATCHED) { 2424 if (_delta2>state(vb) != _delta2>IN_HEAP) { 2425 _delta2>push(vb, _blossom_set>classPrio(vb)  2426 (*_blossom_data)[vb].offset); 2427 } else if ((*_delta2)[vb] > _blossom_set>classPrio(vb)  2428 (*_blossom_data)[vb].offset) { 2429 _delta2>decrease(vb, _blossom_set>classPrio(vb)  2430 (*_blossom_data)[vb].offset); 2431 } 2432 } 2433 } 2434 } 2435 } 2436 } 2437 (*_blossom_data)[blossom].offset = 0; 2438 } 2439 2440 void alternatePath(int even, int tree) { 2441 int odd; 2442 2443 evenToMatched(even, tree); 2444 (*_blossom_data)[even].status = MATCHED; 2445 2446 while ((*_blossom_data)[even].pred != INVALID) { 2447 odd = _blossom_set>find(_ugraph.target((*_blossom_data)[even].pred)); 2448 (*_blossom_data)[odd].status = MATCHED; 2449 oddToMatched(odd); 2450 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; 2451 2452 even = _blossom_set>find(_ugraph.target((*_blossom_data)[odd].pred)); 2453 (*_blossom_data)[even].status = MATCHED; 2454 evenToMatched(even, tree); 2455 (*_blossom_data)[even].next = 2456 _ugraph.oppositeEdge((*_blossom_data)[odd].pred); 2457 } 2458 2459 } 2460 2461 void destroyTree(int tree) { 2462 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { 2463 if ((*_blossom_data)[b].status == EVEN) { 2464 (*_blossom_data)[b].status = MATCHED; 2465 evenToMatched(b, tree); 2466 } else if ((*_blossom_data)[b].status == ODD) { 2467 (*_blossom_data)[b].status = MATCHED; 2468 oddToMatched(b); 2469 } 2470 } 2471 _tree_set>eraseClass(tree); 2472 } 2473 2474 void augmentOnEdge(const UEdge& edge) { 2475 2476 int left = _blossom_set>find(_ugraph.source(edge)); 2477 int right = _blossom_set>find(_ugraph.target(edge)); 2478 2479 int left_tree = _tree_set>find(left); 2480 alternatePath(left, left_tree); 2481 destroyTree(left_tree); 2482 2483 int right_tree = _tree_set>find(right); 2484 alternatePath(right, right_tree); 2485 destroyTree(right_tree); 2486 2487 (*_blossom_data)[left].next = _ugraph.direct(edge, true); 2488 (*_blossom_data)[right].next = _ugraph.direct(edge, false); 2489 } 2490 2491 void extendOnEdge(const Edge& edge) { 2492 int base = _blossom_set>find(_ugraph.target(edge)); 2493 int tree = _tree_set>find(base); 2494 2495 int odd = _blossom_set>find(_ugraph.source(edge)); 2496 _tree_set>insert(odd, tree); 2497 (*_blossom_data)[odd].status = ODD; 2498 matchedToOdd(odd); 2499 (*_blossom_data)[odd].pred = edge; 2500 2501 int even = _blossom_set>find(_ugraph.target((*_blossom_data)[odd].next)); 2502 (*_blossom_data)[even].pred = (*_blossom_data)[even].next; 2503 _tree_set>insert(even, tree); 2504 (*_blossom_data)[even].status = EVEN; 2505 matchedToEven(even, tree); 2506 } 2507 2508 void shrinkOnEdge(const UEdge& uedge, int tree) { 2509 int nca = 1; 2510 std::vector<int> left_path, right_path; 2511 2512 { 2513 std::set<int> left_set, right_set; 2514 int left = _blossom_set>find(_ugraph.source(uedge)); 2515 left_path.push_back(left); 2516 left_set.insert(left); 2517 2518 int right = _blossom_set>find(_ugraph.target(uedge)); 2519 right_path.push_back(right); 2520 right_set.insert(right); 2521 2522 while (true) { 2523 2524 if ((*_blossom_data)[left].pred == INVALID) break; 2525 2526 left = 2527 _blossom_set>find(_ugraph.target((*_blossom_data)[left].pred)); 2528 left_path.push_back(left); 2529 left = 2530 _blossom_set>find(_ugraph.target((*_blossom_data)[left].pred)); 2531 left_path.push_back(left); 2532 2533 left_set.insert(left); 2534 2535 if (right_set.find(left) != right_set.end()) { 2536 nca = left; 2537 break; 2538 } 2539 2540 if ((*_blossom_data)[right].pred == INVALID) break; 2541 2542 right = 2543 _blossom_set>find(_ugraph.target((*_blossom_data)[right].pred)); 2544 right_path.push_back(right); 2545 right = 2546 _blossom_set>find(_ugraph.target((*_blossom_data)[right].pred)); 2547 right_path.push_back(right); 2548 2549 right_set.insert(right); 2550 2551 if (left_set.find(right) != left_set.end()) { 2552 nca = right; 2553 break; 2554 } 2555 2556 } 2557 2558 if (nca == 1) { 2559 if ((*_blossom_data)[left].pred == INVALID) { 2560 nca = right; 2561 while (left_set.find(nca) == left_set.end()) { 2562 nca = 2563 _blossom_set>find(_ugraph.target((*_blossom_data)[nca].pred)); 2564 right_path.push_back(nca); 2565 nca = 2566 _blossom_set>find(_ugraph.target((*_blossom_data)[nca].pred)); 2567 right_path.push_back(nca); 2568 } 2569 } else { 2570 nca = left; 2571 while (right_set.find(nca) == right_set.end()) { 2572 nca = 2573 _blossom_set>find(_ugraph.target((*_blossom_data)[nca].pred)); 2574 left_path.push_back(nca); 2575 nca = 2576 _blossom_set>find(_ugraph.target((*_blossom_data)[nca].pred)); 2577 left_path.push_back(nca); 2578 } 2579 } 2580 } 2581 } 2582 2583 std::vector<int> subblossoms; 2584 Edge prev; 2585 2586 prev = _ugraph.direct(uedge, true); 2587 for (int i = 0; left_path[i] != nca; i += 2) { 2588 subblossoms.push_back(left_path[i]); 2589 (*_blossom_data)[left_path[i]].next = prev; 2590 _tree_set>erase(left_path[i]); 2591 2592 subblossoms.push_back(left_path[i + 1]); 2593 (*_blossom_data)[left_path[i + 1]].status = EVEN; 2594 oddToEven(left_path[i + 1], tree); 2595 _tree_set>erase(left_path[i + 1]); 2596 prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred); 2597 } 2598 2599 int k = 0; 2600 while (right_path[k] != nca) ++k; 2601 2602 subblossoms.push_back(nca); 2603 (*_blossom_data)[nca].next = prev; 2604 2605 for (int i = k  2; i >= 0; i = 2) { 2606 subblossoms.push_back(right_path[i + 1]); 2607 (*_blossom_data)[right_path[i + 1]].status = EVEN; 2608 oddToEven(right_path[i + 1], tree); 2609 _tree_set>erase(right_path[i + 1]); 2610 2611 (*_blossom_data)[right_path[i + 1]].next = 2612 (*_blossom_data)[right_path[i + 1]].pred; 2613 2614 subblossoms.push_back(right_path[i]); 2615 _tree_set>erase(right_path[i]); 2616 } 2617 2618 int surface = 2619 _blossom_set>join(subblossoms.begin(), subblossoms.end()); 2620 2621 for (int i = 0; i < int(subblossoms.size()); ++i) { 2622 if (!_blossom_set>trivial(subblossoms[i])) { 2623 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; 2624 } 2625 (*_blossom_data)[subblossoms[i]].status = MATCHED; 2626 } 2627 2628 (*_blossom_data)[surface].pot = 2 * _delta_sum; 2629 (*_blossom_data)[surface].offset = 0; 2630 (*_blossom_data)[surface].status = EVEN; 2631 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; 2632 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; 2633 2634 _tree_set>insert(surface, tree); 2635 _tree_set>erase(nca); 2636 } 2637 2638 void splitBlossom(int blossom) { 2639 Edge next = (*_blossom_data)[blossom].next; 2640 Edge pred = (*_blossom_data)[blossom].pred; 2641 2642 int tree = _tree_set>find(blossom); 2643 2644 (*_blossom_data)[blossom].status = MATCHED; 2645 oddToMatched(blossom); 2646 if (_delta2>state(blossom) == _delta2>IN_HEAP) { 2647 _delta2>erase(blossom); 2648 } 2649 2650 std::vector<int> subblossoms; 2651 _blossom_set>split(blossom, std::back_inserter(subblossoms)); 2652 2653 Value offset = (*_blossom_data)[blossom].offset; 2654 int b = _blossom_set>find(_ugraph.source(pred)); 2655 int d = _blossom_set>find(_ugraph.source(next)); 2656 2657 int ib, id; 2658 for (int i = 0; i < int(subblossoms.size()); ++i) { 2659 if (subblossoms[i] == b) ib = i; 2660 if (subblossoms[i] == d) id = i; 2661 2662 (*_blossom_data)[subblossoms[i]].offset = offset; 2663 if (!_blossom_set>trivial(subblossoms[i])) { 2664 (*_blossom_data)[subblossoms[i]].pot = 2 * offset; 2665 } 2666 if (_blossom_set>classPrio(subblossoms[i]) != 2667 std::numeric_limits<Value>::max()) { 2668 _delta2>push(subblossoms[i], 2669 _blossom_set>classPrio(subblossoms[i])  2670 (*_blossom_data)[subblossoms[i]].offset); 2671 } 2672 } 2673 2674 if (id > ib ? ((id  ib) % 2 == 0) : ((ib  id) % 2 == 1)) { 2675 for (int i = (id + 1) % subblossoms.size(); 2676 i != ib; i = (i + 2) % subblossoms.size()) { 2677 int sb = subblossoms[i]; 2678 int tb = subblossoms[(i + 1) % subblossoms.size()]; 2679 (*_blossom_data)[sb].next = 2680 _ugraph.oppositeEdge((*_blossom_data)[tb].next); 2681 } 2682 2683 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { 2684 int sb = subblossoms[i]; 2685 int tb = subblossoms[(i + 1) % subblossoms.size()]; 2686 int ub = subblossoms[(i + 2) % subblossoms.size()]; 2687 2688 (*_blossom_data)[sb].status = ODD; 2689 matchedToOdd(sb); 2690 _tree_set>insert(sb, tree); 2691 (*_blossom_data)[sb].pred = pred; 2692 (*_blossom_data)[sb].next = 2693 _ugraph.oppositeEdge((*_blossom_data)[tb].next); 2694 2695 pred = (*_blossom_data)[ub].next; 2696 2697 (*_blossom_data)[tb].status = EVEN; 2698 matchedToEven(tb, tree); 2699 _tree_set>insert(tb, tree); 2700 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; 2701 } 2702 2703 (*_blossom_data)[subblossoms[id]].status = ODD; 2704 matchedToOdd(subblossoms[id]); 2705 _tree_set>insert(subblossoms[id], tree); 2706 (*_blossom_data)[subblossoms[id]].next = next; 2707 (*_blossom_data)[subblossoms[id]].pred = pred; 2708 2709 } else { 2710 2711 for (int i = (ib + 1) % subblossoms.size(); 2712 i != id; i = (i + 2) % subblossoms.size()) { 2713 int sb = subblossoms[i]; 2714 int tb = subblossoms[(i + 1) % subblossoms.size()]; 2715 (*_blossom_data)[sb].next = 2716 _ugraph.oppositeEdge((*_blossom_data)[tb].next); 2717 } 2718 2719 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { 2720 int sb = subblossoms[i]; 2721 int tb = subblossoms[(i + 1) % subblossoms.size()]; 2722 int ub = subblossoms[(i + 2) % subblossoms.size()]; 2723 2724 (*_blossom_data)[sb].status = ODD; 2725 matchedToOdd(sb); 2726 _tree_set>insert(sb, tree); 2727 (*_blossom_data)[sb].next = next; 2728 (*_blossom_data)[sb].pred = 2729 _ugraph.oppositeEdge((*_blossom_data)[tb].next); 2730 2731 (*_blossom_data)[tb].status = EVEN; 2732 matchedToEven(tb, tree); 2733 _tree_set>insert(tb, tree); 2734 (*_blossom_data)[tb].pred = 2735 (*_blossom_data)[tb].next = 2736 _ugraph.oppositeEdge((*_blossom_data)[ub].next); 2737 next = (*_blossom_data)[ub].next; 2738 } 2739 2740 (*_blossom_data)[subblossoms[ib]].status = ODD; 2741 matchedToOdd(subblossoms[ib]); 2742 _tree_set>insert(subblossoms[ib], tree); 2743 (*_blossom_data)[subblossoms[ib]].next = next; 2744 (*_blossom_data)[subblossoms[ib]].pred = pred; 2745 } 2746 _tree_set>erase(blossom); 2747 } 2748 2749 void extractBlossom(int blossom, const Node& base, const Edge& matching) { 2750 if (_blossom_set>trivial(blossom)) { 2751 int bi = (*_node_index)[base]; 2752 Value pot = (*_node_data)[bi].pot; 2753 2754 _matching>set(base, matching); 2755 _blossom_node_list.push_back(base); 2756 _node_potential>set(base, pot); 2757 } else { 2758 2759 Value pot = (*_blossom_data)[blossom].pot; 2760 int bn = _blossom_node_list.size(); 2761 2762 std::vector<int> subblossoms; 2763 _blossom_set>split(blossom, std::back_inserter(subblossoms)); 2764 int b = _blossom_set>find(base); 2765 int ib = 1; 2766 for (int i = 0; i < int(subblossoms.size()); ++i) { 2767 if (subblossoms[i] == b) { ib = i; break; } 2768 } 2769 2770 for (int i = 1; i < int(subblossoms.size()); i += 2) { 2771 int sb = subblossoms[(ib + i) % subblossoms.size()]; 2772 int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; 2773 2774 Edge m = (*_blossom_data)[tb].next; 2775 extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m)); 2776 extractBlossom(tb, _ugraph.source(m), m); 2777 } 2778 extractBlossom(subblossoms[ib], base, matching); 2779 2780 int en = _blossom_node_list.size(); 2781 2782 _blossom_potential.push_back(BlossomVariable(bn, en, pot)); 2783 } 2784 } 2785 2786 void extractMatching() { 2787 std::vector<int> blossoms; 2788 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { 2789 blossoms.push_back(c); 2790 } 2791 2792 for (int i = 0; i < int(blossoms.size()); ++i) { 2793 2794 Value offset = (*_blossom_data)[blossoms[i]].offset; 2795 (*_blossom_data)[blossoms[i]].pot += 2 * offset; 2796 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); 2797 n != INVALID; ++n) { 2798 (*_node_data)[(*_node_index)[n]].pot = offset; 2799 } 2800 2801 Edge matching = (*_blossom_data)[blossoms[i]].next; 2802 Node base = _ugraph.source(matching); 2803 extractBlossom(blossoms[i], base, matching); 2804 } 2805 } 2806 2807 public: 2808 2809 /// \brief Constructor 2810 /// 2811 /// Constructor. 2812 MaxWeightedPerfectMatching(const UGraph& ugraph, const WeightMap& weight) 2813 : _ugraph(ugraph), _weight(weight), _matching(0), 2814 _node_potential(0), _blossom_potential(), _blossom_node_list(), 2815 _node_num(0), _blossom_num(0), 2816 2817 _blossom_index(0), _blossom_set(0), _blossom_data(0), 2818 _node_index(0), _node_heap_index(0), _node_data(0), 2819 _tree_set_index(0), _tree_set(0), 2820 2821 _delta2_index(0), _delta2(0), 2822 _delta3_index(0), _delta3(0), 2823 _delta4_index(0), _delta4(0), 2824 2825 _delta_sum() {} 2826 2827 ~MaxWeightedPerfectMatching() { 2828 destroyStructures(); 2829 } 2830 2831 /// \name Execution control 2832 /// The simplest way to execute the algorithm is to use the member 2833 /// \c run() member function. 2834 2835 ///@{ 2836 2837 /// \brief Initialize the algorithm 2838 /// 2839 /// Initialize the algorithm 2840 void init() { 2841 createStructures(); 2842 2843 for (EdgeIt e(_ugraph); e != INVALID; ++e) { 2844 _node_heap_index>set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP); 2845 } 2846 for (UEdgeIt e(_ugraph); e != INVALID; ++e) { 2847 _delta3_index>set(e, _delta3>PRE_HEAP); 2848 } 2849 for (int i = 0; i < _blossom_num; ++i) { 2850 _delta2_index>set(i, _delta2>PRE_HEAP); 2851 _delta4_index>set(i, _delta4>PRE_HEAP); 2852 } 2853 2854 int index = 0; 2855 for (NodeIt n(_ugraph); n != INVALID; ++n) { 2856 Value max = std::numeric_limits<Value>::min(); 2857 for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) { 2858 if (_ugraph.target(e) == n) continue; 2859 if ((dualScale * _weight[e]) / 2 > max) { 2860 max = (dualScale * _weight[e]) / 2; 2861 } 2862 } 2863 _node_index>set(n, index); 2864 (*_node_data)[index].pot = max; 2865 int blossom = 2866 _blossom_set>insert(n, std::numeric_limits<Value>::max()); 2867 2868 _tree_set>insert(blossom); 2869 2870 (*_blossom_data)[blossom].status = EVEN; 2871 (*_blossom_data)[blossom].pred = INVALID; 2872 (*_blossom_data)[blossom].next = INVALID; 2873 (*_blossom_data)[blossom].pot = 0; 2874 (*_blossom_data)[blossom].offset = 0; 2875 ++index; 2876 } 2877 for (UEdgeIt e(_ugraph); e != INVALID; ++e) { 2878 int si = (*_node_index)[_ugraph.source(e)]; 2879 int ti = (*_node_index)[_ugraph.target(e)]; 2880 if (_ugraph.source(e) != _ugraph.target(e)) { 2881 _delta3>push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot  2882 dualScale * _weight[e]) / 2); 2883 } 2884 } 2885 } 2886 2887 /// \brief Starts the algorithm 2888 /// 2889 /// Starts the algorithm 2890 bool start() { 2891 enum OpType { 2892 D2, D3, D4 2893 }; 2894 2895 int unmatched = _node_num; 2896 while (unmatched > 0) { 2897 Value d2 = !_delta2>empty() ? 2898 _delta2>prio() : std::numeric_limits<Value>::max(); 2899 2900 Value d3 = !_delta3>empty() ? 2901 _delta3>prio() : std::numeric_limits<Value>::max(); 2902 2903 Value d4 = !_delta4>empty() ? 2904 _delta4>prio() : std::numeric_limits<Value>::max(); 2905 2906 _delta_sum = d2; OpType ot = D2; 2907 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } 2908 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } 2909 2910 if (_delta_sum == std::numeric_limits<Value>::max()) { 2911 return false; 2912 } 2913 2914 switch (ot) { 2915 case D2: 2916 { 2917 int blossom = _delta2>top(); 2918 Node n = _blossom_set>classTop(blossom); 2919 Edge e = (*_node_data)[(*_node_index)[n]].heap.top(); 2920 extendOnEdge(e); 2921 } 2922 break; 2923 case D3: 2924 { 2925 UEdge e = _delta3>top(); 2926 2927 int left_blossom = _blossom_set>find(_ugraph.source(e)); 2928 int right_blossom = _blossom_set>find(_ugraph.target(e)); 2929 2930 if (left_blossom == right_blossom) { 2931 _delta3>pop(); 2932 } else { 2933 int left_tree = _tree_set>find(left_blossom); 2934 int right_tree = _tree_set>find(right_blossom); 2935 2936 if (left_tree == right_tree) { 2937 shrinkOnEdge(e, left_tree); 2938 } else { 2939 augmentOnEdge(e); 2940 unmatched = 2; 2941 } 2942 } 2943 } break; 2944 case D4: 2945 splitBlossom(_delta4>top()); 2946 break; 2947 } 2948 } 2949 extractMatching(); 2950 return true; 2951 } 2952 2953 /// \brief Runs %MaxWeightedPerfectMatching algorithm. 2954 /// 2955 /// This method runs the %MaxWeightedPerfectMatching algorithm. 2956 /// 2957 /// \note mwm.run() is just a shortcut of the following code. 2958 /// \code 2959 /// mwm.init(); 2960 /// mwm.start(); 2961 /// \endcode 2962 bool run() { 2963 init(); 2964 return start(); 2965 } 2966 2967 /// @} 2968 2969 /// \name Primal solution 2970 /// Functions for get the primal solution, ie. the matching. 2971 2972 /// @{ 2973 2974 /// \brief Returns the matching value. 2975 /// 2976 /// Returns the matching value. 2977 Value matchingValue() const { 2978 Value sum = 0; 2979 for (NodeIt n(_ugraph); n != INVALID; ++n) { 2980 if ((*_matching)[n] != INVALID) { 2981 sum += _weight[(*_matching)[n]]; 2982 } 2983 } 2984 return sum /= 2; 2985 } 2986 2987 /// \brief Returns true when the edge is in the matching. 2988 /// 2989 /// Returns true when the edge is in the matching. 2990 bool matching(const UEdge& edge) const { 2991 return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true); 2992 } 2993 2994 /// \brief Returns the incident matching edge. 2995 /// 2996 /// Returns the incident matching edge from given node. 2997 Edge matching(const Node& node) const { 2998 return (*_matching)[node]; 2999 } 3000 3001 /// \brief Returns the mate of the node. 3002 /// 3003 /// Returns the adjancent node in a mathcing edge. 3004 Node mate(const Node& node) const { 3005 return _ugraph.target((*_matching)[node]); 3006 } 3007 3008 /// @} 3009 3010 /// \name Dual solution 3011 /// Functions for get the dual solution. 3012 3013 /// @{ 3014 3015 /// \brief Returns the value of the dual solution. 3016 /// 3017 /// Returns the value of the dual solution. It should be equal to 3018 /// the primal value scaled by \ref dualScale "dual scale". 3019 Value dualValue() const { 3020 Value sum = 0; 3021 for (NodeIt n(_ugraph); n != INVALID; ++n) { 3022 sum += nodeValue(n); 3023 } 3024 for (int i = 0; i < blossomNum(); ++i) { 3025 sum += blossomValue(i) * (blossomSize(i) / 2); 3026 } 3027 return sum; 3028 } 3029 3030 /// \brief Returns the value of the node. 3031 /// 3032 /// Returns the the value of the node. 3033 Value nodeValue(const Node& n) const { 3034 return (*_node_potential)[n]; 3035 } 3036 3037 /// \brief Returns the number of the blossoms in the basis. 3038 /// 3039 /// Returns the number of the blossoms in the basis. 3040 /// \see BlossomIt 3041 int blossomNum() const { 3042 return _blossom_potential.size(); 3043 } 3044 3045 3046 /// \brief Returns the number of the nodes in the blossom. 3047 /// 3048 /// Returns the number of the nodes in the blossom. 3049 int blossomSize(int k) const { 3050 return _blossom_potential[k].end  _blossom_potential[k].begin; 3051 } 3052 3053 /// \brief Returns the value of the blossom. 3054 /// 3055 /// Returns the the value of the blossom. 3056 /// \see BlossomIt 3057 Value blossomValue(int k) const { 3058 return _blossom_potential[k].value; 3059 } 3060 3061 /// \brief Lemon iterator for get the items of the blossom. 3062 /// 3063 /// Lemon iterator for get the nodes of the blossom. This class 3064 /// provides a common style lemon iterator which gives back a 3065 /// subset of the nodes. 3066 class BlossomIt { 3067 public: 3068 3069 /// \brief Constructor. 3070 /// 3071 /// Constructor for get the nodes of the variable. 3072 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) 3073 : _algorithm(&algorithm) 3074 { 3075 _index = _algorithm>_blossom_potential[variable].begin; 3076 _last = _algorithm>_blossom_potential[variable].end; 3077 } 3078 3079 /// \brief Invalid constructor. 3080 /// 3081 /// Invalid constructor. 3082 BlossomIt(Invalid) : _index(1) {} 3083 3084 /// \brief Conversion to node. 3085 /// 3086 /// Conversion to node. 3087 operator Node() const { 3088 return _algorithm ? _algorithm>_blossom_node_list[_index] : INVALID; 3089 } 3090 3091 /// \brief Increment operator. 3092 /// 3093 /// Increment operator. 3094 BlossomIt& operator++() { 3095 ++_index; 3096 if (_index == _last) { 3097 _index = 1; 3098 } 3099 return *this; 3100 } 3101 3102 bool operator==(const BlossomIt& it) const { 3103 return _index == it._index; 3104 } 3105 bool operator!=(const BlossomIt& it) const { 3106 return _index != it._index; 3107 } 3108 3109 private: 3110 const MaxWeightedPerfectMatching* _algorithm; 3111 int _last; 3112 int _index; 3113 }; 3114 3115 /// @} 3116 3117 }; 3118 621 3119 622 3120 } //END OF NAMESPACE LEMON 
lemon/unionfind.h
r2542 r2548 925 925 }; 926 926 927 /// \ingroup auxdat 928 /// 929 /// \brief A \e UnionFind data structure implementation which 930 /// is able to store a priority for each item and retrieve the minimum of 931 /// each class. 932 /// 933 /// A \e UnionFind data structure implementation which is able to 934 /// store a priority for each item and retrieve the minimum of each 935 /// class. In addition, it supports the joining and splitting the 936 /// components. If you don't need this feature then you makes 937 /// better to use the \ref UnionFind class which is more efficient. 938 /// 939 /// The unionfind data strcuture based on a (2, 16)tree with a 940 /// tournament minimum selection on the internal nodes. The insert 941 /// operation takes O(1), the find, set, decrease and increase takes 942 /// O(log(n)), where n is the number of nodes in the current 943 /// component. The complexity of join and split is O(log(n)*k), 944 /// where n is the sum of the number of the nodes and k is the 945 /// number of joined components or the number of the components 946 /// after the split. 947 /// 948 /// \pre You need to add all the elements by the \ref insert() 949 /// method. 950 /// 951 template <typename _Value, typename _ItemIntMap, 952 typename _Comp = std::less<_Value> > 953 class HeapUnionFind { 954 public: 955 956 typedef _Value Value; 957 typedef typename _ItemIntMap::Key Item; 958 959 typedef _ItemIntMap ItemIntMap; 960 961 typedef _Comp Comp; 962 963 private: 964 965 static const int cmax = 3; 966 967 ItemIntMap& index; 968 969 struct ClassNode { 970 int parent; 971 int depth; 972 973 int left, right; 974 int next, prev; 975 }; 976 977 int first_class; 978 int first_free_class; 979 std::vector<ClassNode> classes; 980 981 int newClass() { 982 if (first_free_class < 0) { 983 int id = classes.size(); 984 classes.push_back(ClassNode()); 985 return id; 986 } else { 987 int id = first_free_class; 988 first_free_class = classes[id].next; 989 return id; 990 } 991 } 992 993 void deleteClass(int id) { 994 classes[id].next = first_free_class; 995 first_free_class = id; 996 } 997 998 struct ItemNode { 999 int parent; 1000 Item item; 1001 Value prio; 1002 int next, prev; 1003 int left, right; 1004 int size; 1005 }; 1006 1007 int first_free_node; 1008 std::vector<ItemNode> nodes; 1009 1010 int newNode() { 1011 if (first_free_node < 0) { 1012 int id = nodes.size(); 1013 nodes.push_back(ItemNode()); 1014 return id; 1015 } else { 1016 int id = first_free_node; 1017 first_free_node = nodes[id].next; 1018 return id; 1019 } 1020 } 1021 1022 void deleteNode(int id) { 1023 nodes[id].next = first_free_node; 1024 first_free_node = id; 1025 } 1026 1027 Comp comp; 1028 1029 int findClass(int id) const { 1030 int kd = id; 1031 while (kd >= 0) { 1032 kd = nodes[kd].parent; 1033 } 1034 return ~kd; 1035 } 1036 1037 int leftNode(int id) const { 1038 int kd = ~(classes[id].parent); 1039 for (int i = 0; i < classes[id].depth; ++i) { 1040 kd = nodes[kd].left; 1041 } 1042 return kd; 1043 } 1044 1045 int nextNode(int id) const { 1046 int depth = 0; 1047 while (id >= 0 && nodes[id].next == 1) { 1048 id = nodes[id].parent; 1049 ++depth; 1050 } 1051 if (id < 0) { 1052 return 1; 1053 } 1054 id = nodes[id].next; 1055 while (depth) { 1056 id = nodes[id].left; 1057 } 1058 return id; 1059 } 1060 1061 1062 void setPrio(int id) { 1063 int jd = nodes[id].left; 1064 nodes[id].prio = nodes[jd].prio; 1065 nodes[id].item = nodes[jd].item; 1066 jd = nodes[jd].next; 1067 while (jd != 1) { 1068 if (comp(nodes[jd].prio, nodes[id].prio)) { 1069 nodes[id].prio = nodes[jd].prio; 1070 nodes[id].item = nodes[jd].item; 1071 } 1072 jd = nodes[jd].next; 1073 } 1074 } 1075 1076 void push(int id, int jd) { 1077 nodes[id].size = 1; 1078 nodes[id].left = nodes[id].right = jd; 1079 nodes[jd].next = nodes[jd].prev = 1; 1080 nodes[jd].parent = id; 1081 } 1082 1083 void pushAfter(int id, int jd) { 1084 int kd = nodes[id].parent; 1085 if (nodes[id].next != 1) { 1086 nodes[nodes[id].next].prev = jd; 1087 if (kd >= 0) { 1088 nodes[kd].size += 1; 1089 } 1090 } else { 1091 if (kd >= 0) { 1092 nodes[kd].right = jd; 1093 nodes[kd].size += 1; 1094 } 1095 } 1096 nodes[jd].next = nodes[id].next; 1097 nodes[jd].prev = id; 1098 nodes[id].next = jd; 1099 nodes[jd].parent = kd; 1100 } 1101 1102 void pushRight(int id, int jd) { 1103 nodes[id].size += 1; 1104 nodes[jd].prev = nodes[id].right; 1105 nodes[jd].next = 1; 1106 nodes[nodes[id].right].next = jd; 1107 nodes[id].right = jd; 1108 nodes[jd].parent = id; 1109 } 1110 1111 void popRight(int id) { 1112 nodes[id].size = 1; 1113 int jd = nodes[id].right; 1114 nodes[nodes[jd].prev].next = 1; 1115 nodes[id].right = nodes[jd].prev; 1116 } 1117 1118 void splice(int id, int jd) { 1119 nodes[id].size += nodes[jd].size; 1120 nodes[nodes[id].right].next = nodes[jd].left; 1121 nodes[nodes[jd].left].prev = nodes[id].right; 1122 int kd = nodes[jd].left; 1123 while (kd != 1) { 1124 nodes[kd].parent = id; 1125 kd = nodes[kd].next; 1126 } 1127 } 1128 1129 void split(int id, int jd) { 1130 int kd = nodes[id].parent; 1131 nodes[kd].right = nodes[id].prev; 1132 nodes[nodes[id].prev].next = 1; 1133 1134 nodes[jd].left = id; 1135 nodes[id].prev = 1; 1136 int num = 0; 1137 while (id != 1) { 1138 nodes[id].parent = jd; 1139 nodes[jd].right = id; 1140 id = nodes[id].next; 1141 ++num; 1142 } 1143 nodes[kd].size = num; 1144 nodes[jd].size = num; 1145 } 1146 1147 void pushLeft(int id, int jd) { 1148 nodes[id].size += 1; 1149 nodes[jd].next = nodes[id].left; 1150 nodes[jd].prev = 1; 1151 nodes[nodes[id].left].prev = jd; 1152 nodes[id].left = jd; 1153 nodes[jd].parent = id; 1154 } 1155 1156 void popLeft(int id) { 1157 nodes[id].size = 1; 1158 int jd = nodes[id].left; 1159 nodes[nodes[jd].next].prev = 1; 1160 nodes[id].left = nodes[jd].next; 1161 } 1162 1163 void repairLeft(int id) { 1164 int jd = ~(classes[id].parent); 1165 while (nodes[jd].left != 1) { 1166 int kd = nodes[jd].left; 1167 if (nodes[jd].size == 1) { 1168 if (nodes[jd].parent < 0) { 1169 classes[id].parent = ~kd; 1170 classes[id].depth = 1; 1171 nodes[kd].parent = ~id; 1172 deleteNode(jd); 1173 jd = kd; 1174 } else { 1175 int pd = nodes[jd].parent; 1176 if (nodes[nodes[jd].next].size < cmax) { 1177 pushLeft(nodes[jd].next, nodes[jd].left); 1178 if (less(nodes[jd].left, nodes[jd].next)) { 1179 nodes[nodes[jd].next].prio = nodes[nodes[jd].left].prio; 1180 nodes[nodes[jd].next].item = nodes[nodes[jd].left].item; 1181 } 1182 popLeft(pd); 1183 deleteNode(jd); 1184 jd = pd; 1185 } else { 1186 int ld = nodes[nodes[jd].next].left; 1187 popLeft(nodes[jd].next); 1188 pushRight(jd, ld); 1189 if (less(ld, nodes[jd].left)) { 1190 nodes[jd].item = nodes[ld].item; 1191 nodes[jd].prio = nodes[jd].prio; 1192 } 1193 if (nodes[nodes[jd].next].item == nodes[ld].item) { 1194 setPrio(nodes[jd].next); 1195 } 1196 jd = nodes[jd].left; 1197 } 1198 } 1199 } else { 1200 jd = nodes[jd].left; 1201 } 1202 } 1203 } 1204 1205 void repairRight(int id) { 1206 int jd = ~(classes[id].parent); 1207 while (nodes[jd].right != 1) { 1208 int kd = nodes[jd].right; 1209 if (nodes[jd].size == 1) { 1210 if (nodes[jd].parent < 0) { 1211 classes[id].parent = ~kd; 1212 classes[id].depth = 1; 1213 nodes[kd].parent = ~id; 1214 deleteNode(jd); 1215 jd = kd; 1216 } else { 1217 int pd = nodes[jd].parent; 1218 if (nodes[nodes[jd].prev].size < cmax) { 1219 pushRight(nodes[jd].prev, nodes[jd].right); 1220 if (less(nodes[jd].right, nodes[jd].prev)) { 1221 nodes[nodes[jd].prev].prio = nodes[nodes[jd].right].prio; 1222 nodes[nodes[jd].prev].item = nodes[nodes[jd].right].item; 1223 } 1224 popRight(pd); 1225 deleteNode(jd); 1226 jd = pd; 1227 } else { 1228 int ld = nodes[nodes[jd].prev].right; 1229 popRight(nodes[jd].prev); 1230 pushLeft(jd, ld); 1231 if (less(ld, nodes[jd].right)) { 1232 nodes[jd].item = nodes[ld].item; 1233 nodes[jd].prio = nodes[jd].prio; 1234 } 1235 if (nodes[nodes[jd].prev].item == nodes[ld].item) { 1236 setPrio(nodes[jd].prev); 1237 } 1238 jd = nodes[jd].right; 1239 } 1240 } 1241 } else { 1242 jd = nodes[jd].right; 1243 } 1244 } 1245 } 1246 1247 1248 bool less(int id, int jd) const { 1249 return comp(nodes[id].prio, nodes[jd].prio); 1250 } 1251 1252 bool equal(int id, int jd) const { 1253 return !less(id, jd) && !less(jd, id); 1254 } 1255 1256 1257 public: 1258 1259 /// \brief Returns true when the given class is alive. 1260 /// 1261 /// Returns true when the given class is alive, ie. the class is 1262 /// not nested into other class. 1263 bool alive(int cls) const { 1264 return classes[cls].parent < 0; 1265 } 1266 1267 /// \brief Returns true when the given class is trivial. 1268 /// 1269 /// Returns true when the given class is trivial, ie. the class 1270 /// contains just one item directly. 1271 bool trivial(int cls) const { 1272 return classes[cls].left == 1; 1273 } 1274 1275 /// \brief Constructs the unionfind. 1276 /// 1277 /// Constructs the unionfind. 1278 /// \brief _index The index map of the unionfind. The data 1279 /// structure uses internally for store references. 1280 HeapUnionFind(ItemIntMap& _index) 1281 : index(_index), first_class(1), 1282 first_free_class(1), first_free_node(1) {} 1283 1284 /// \brief Insert a new node into a new component. 1285 /// 1286 /// Insert a new node into a new component. 1287 /// \param item The item of the new node. 1288 /// \param prio The priority of the new node. 1289 /// \return The class id of the oneitemheap. 1290 int insert(const Item& item, const Value& prio) { 1291 int id = newNode(); 1292 nodes[id].item = item; 1293 nodes[id].prio = prio; 1294 nodes[id].size = 0; 1295 1296 nodes[id].prev = 1; 1297 nodes[id].next = 1; 1298 1299 nodes[id].left = 1; 1300 nodes[id].right = 1; 1301 1302 nodes[id].item = item; 1303 index[item] = id; 1304 1305 int class_id = newClass(); 1306 classes[class_id].parent = ~id; 1307 classes[class_id].depth = 0; 1308 1309 classes[class_id].left = 1; 1310 classes[class_id].right = 1; 1311 1312 if (first_class != 1) { 1313 classes[first_class].prev = class_id; 1314 } 1315 classes[class_id].next = first_class; 1316 classes[class_id].prev = 1; 1317 first_class = class_id; 1318 1319 nodes[id].parent = ~class_id; 1320 1321 return class_id; 1322 } 1323 1324 /// \brief The class of the item. 1325 /// 1326 /// \return The alive class id of the item, which is not nested into 1327 /// other classes. 1328 /// 1329 /// The time complexity is O(log(n)). 1330 int find(const Item& item) const { 1331 return findClass(index[item]); 1332 } 1333 1334 /// \brief Joins the classes. 1335 /// 1336 /// The current function joins the given classes. The parameter is 1337 /// an STL range which should be contains valid class ids. The 1338 /// time complexity is O(log(n)*k) where n is the overall number 1339 /// of the joined nodes and k is the number of classes. 1340 /// \return The class of the joined classes. 1341 /// \pre The range should contain at least two class ids. 1342 template <typename Iterator> 1343 int join(Iterator begin, Iterator end) { 1344 std::vector<int> cs; 1345 for (Iterator it = begin; it != end; ++it) { 1346 cs.push_back(*it); 1347 } 1348 1349 int class_id = newClass(); 1350 { // creation unionfind 1351 1352 if (first_class != 1) { 1353 classes[first_class].prev = class_id; 1354 } 1355 classes[class_id].next = first_class; 1356 classes[class_id].prev = 1; 1357 first_class = class_id; 1358 1359 classes[class_id].depth = classes[cs[0]].depth; 1360 classes[class_id].parent = classes[cs[0]].parent; 1361 nodes[~(classes[class_id].parent)].parent = ~class_id; 1362 1363 int l = cs[0]; 1364 1365 classes[class_id].left = l; 1366 classes[class_id].right = l; 1367 1368 if (classes[l].next != 1) { 1369 classes[classes[l].next].prev = classes[l].prev; 1370 } 1371 classes[classes[l].prev].next = classes[l].next; 1372 1373 classes[l].prev = 1; 1374 classes[l].next = 1; 1375 1376 classes[l].depth = leftNode(l); 1377 classes[l].parent = class_id; 1378 1379 } 1380 1381 { // merging of heap 1382 int l = class_id; 1383 for (int ci = 1; ci < int(cs.size()); ++ci) { 1384 int r = cs[ci]; 1385 int rln = leftNode(r); 1386 if (classes[l].depth > classes[r].depth) { 1387 int id = ~(classes[l].parent); 1388 for (int i = classes[r].depth + 1; i < classes[l].depth; ++i) { 1389 id = nodes[id].right; 1390 } 1391 while (id >= 0 && nodes[id].size == cmax) { 1392 int new_id = newNode(); 1393 int right_id = nodes[id].right; 1394 1395 popRight(id); 1396 if (nodes[id].item == nodes[right_id].item) { 1397 setPrio(id); 1398 } 1399 push(new_id, right_id); 1400 pushRight(new_id, ~(classes[r].parent)); 1401 setPrio(new_id); 1402 1403 id = nodes[id].parent; 1404 classes[r].parent = ~new_id; 1405 } 1406 if (id < 0) { 1407 int new_parent = newNode(); 1408 nodes[new_parent].next = 1; 1409 nodes[new_parent].prev = 1; 1410 nodes[new_parent].parent = ~l; 1411 1412 push(new_parent, ~(classes[l].parent)); 1413 pushRight(new_parent, ~(classes[r].parent)); 1414 setPrio(new_parent); 1415 1416 classes[l].parent = ~new_parent; 1417 classes[l].depth += 1; 1418 } else { 1419 pushRight(id, ~(classes[r].parent)); 1420 while (id >= 0 && less(~(classes[r].parent), id)) { 1421 nodes[id].prio = nodes[~(classes[r].parent)].prio; 1422 nodes[id].item = nodes[~(classes[r].parent)].item; 1423 id = nodes[id].parent; 1424 } 1425 } 1426 } else if (classes[r].depth > classes[l].depth) { 1427 int id = ~(classes[r].parent); 1428 for (int i = classes[l].depth + 1; i < classes[r].depth; ++i) { 1429 id = nodes[id].left; 1430 } 1431 while (id >= 0 && nodes[id].size == cmax) { 1432 int new_id = newNode(); 1433 int left_id = nodes[id].left; 1434 1435 popLeft(id); 1436 if (nodes[id].prio == nodes[left_id].prio) { 1437 setPrio(id); 1438 } 1439 push(new_id, left_id); 1440 pushLeft(new_id, ~(classes[l].parent)); 1441 setPrio(new_id); 1442 1443 id = nodes[id].parent; 1444 classes[l].parent = ~new_id; 1445 1446 } 1447 if (id < 0) { 1448 int new_parent = newNode(); 1449 nodes[new_parent].next = 1; 1450 nodes[new_parent].prev = 1; 1451 nodes[new_parent].parent = ~l; 1452 1453 push(new_parent, ~(classes[r].parent)); 1454 pushLeft(new_parent, ~(classes[l].parent)); 1455 setPrio(new_parent); 1456 1457 classes[r].parent = ~new_parent; 1458 classes[r].depth += 1; 1459 } else { 1460 pushLeft(id, ~(classes[l].parent)); 1461 while (id >= 0 && less(~(classes[l].parent), id)) { 1462 nodes[id].prio = nodes[~(classes[l].parent)].prio; 1463 nodes[id].item = nodes[~(classes[l].parent)].item; 1464 id = nodes[id].parent; 1465 } 1466 } 1467 nodes[~(classes[r].parent)].parent = ~l; 1468 classes[l].parent = classes[r].parent; 1469 classes[l].depth = classes[r].depth; 1470 } else { 1471 if (classes[l].depth != 0 && 1472 nodes[~(classes[l].parent)].size + 1473 nodes[~(classes[r].parent)].size <= cmax) { 1474 splice(~(classes[l].parent), ~(classes[r].parent)); 1475 deleteNode(~(classes[r].parent)); 1476 if (less(~(classes[r].parent), ~(classes[l].parent))) { 1477 nodes[~(classes[l].parent)].prio = 1478 nodes[~(classes[r].parent)].prio; 1479 nodes[~(classes[l].parent)].item = 1480 nodes[~(classes[r].parent)].item; 1481 } 1482 } else { 1483 int new_parent = newNode(); 1484 nodes[new_parent].next = nodes[new_parent].prev = 1; 1485 push(new_parent, ~(classes[l].parent)); 1486 pushRight(new_parent, ~(classes[r].parent)); 1487 setPrio(new_parent); 1488 1489 classes[l].parent = ~new_parent; 1490 classes[l].depth += 1; 1491 nodes[new_parent].parent = ~l; 1492 } 1493 } 1494 if (classes[r].next != 1) { 1495 classes[classes[r].next].prev = classes[r].prev; 1496 } 1497 classes[classes[r].prev].next = classes[r].next; 1498 1499 classes[r].prev = classes[l].right; 1500 classes[classes[l].right].next = r; 1501 classes[l].right = r; 1502 classes[r].parent = l; 1503 1504 classes[r].next = 1; 1505 classes[r].depth = rln; 1506 } 1507 } 1508 return class_id; 1509 } 1510 1511 /// \brief Split the class to subclasses. 1512 /// 1513 /// The current function splits the given class. The join, which 1514 /// made the current class, stored a reference to the 1515 /// subclasses. The \c splitClass() member restores the classes 1516 /// and creates the heaps. The parameter is an STL output iterator 1517 /// which will be filled with the subclass ids. The time 1518 /// complexity is O(log(n)*k) where n is the overall number of 1519 /// nodes in the splitted classes and k is the number of the 1520 /// classes. 1521 template <typename Iterator> 1522 void split(int cls, Iterator out) { 1523 std::vector<int> cs; 1524 { // splitting unionfind 1525 int id = cls; 1526 int l = classes[id].left; 1527 1528 classes[l].parent = classes[id].parent; 1529 classes[l].depth = classes[id].depth; 1530 1531 nodes[~(classes[l].parent)].parent = ~l; 1532 1533 *out++ = l; 1534 1535 while (l != 1) { 1536 cs.push_back(l); 1537 l = classes[l].next; 1538 } 1539 1540 classes[classes[id].right].next = first_class; 1541 classes[first_class].prev = classes[id].right; 1542 first_class = classes[id].left; 1543 1544 if (classes[id].next != 1) { 1545 classes[classes[id].next].prev = classes[id].prev; 1546 } 1547 classes[classes[id].prev].next = classes[id].next; 1548 1549 deleteClass(id); 1550 } 1551 1552 { 1553 for (int i = 1; i < int(cs.size()); ++i) { 1554 int l = classes[cs[i]].depth; 1555 while (nodes[nodes[l].parent].left == l) { 1556 l = nodes[l].parent; 1557 } 1558 int r = l; 1559 while (nodes[l].parent >= 0) { 1560 l = nodes[l].parent; 1561 int new_node = newNode(); 1562 1563 nodes[new_node].prev = 1; 1564 nodes[new_node].next = 1; 1565 1566 split(r, new_node); 1567 pushAfter(l, new_node); 1568 setPrio(l); 1569 setPrio(new_node); 1570 r = new_node; 1571 } 1572 classes[cs[i]].parent = ~r; 1573 classes[cs[i]].depth = classes[~(nodes[l].parent)].depth; 1574 nodes[r].parent = ~cs[i]; 1575 1576 nodes[l].next = 1; 1577 nodes[r].prev = 1; 1578 1579 repairRight(~(nodes[l].parent)); 1580 repairLeft(cs[i]); 1581 1582 *out++ = cs[i]; 1583 } 1584 } 1585 } 1586 1587 /// \brief Gives back the priority of the current item. 1588 /// 1589 /// \return Gives back the priority of the current item. 1590 const Value& operator[](const Item& item) const { 1591 return nodes[index[item]].prio; 1592 } 1593 1594 /// \brief Sets the priority of the current item. 1595 /// 1596 /// Sets the priority of the current item. 1597 void set(const Item& item, const Value& prio) { 1598 if (comp(prio, nodes[index[item]].prio)) { 1599 decrease(item, prio); 1600 } else if (!comp(prio, nodes[index[item]].prio)) { 1601 increase(item, prio); 1602 } 1603 } 1604 1605 /// \brief Increase the priority of the current item. 1606 /// 1607 /// Increase the priority of the current item. 1608 void increase(const Item& item, const Value& prio) { 1609 int id = index[item]; 1610 int kd = nodes[id].parent; 1611 nodes[id].prio = prio; 1612 while (kd >= 0 && nodes[kd].item == item) { 1613 setPrio(kd); 1614 kd = nodes[kd].parent; 1615 } 1616 } 1617 1618 /// \brief Increase the priority of the current item. 1619 /// 1620 /// Increase the priority of the current item. 1621 void decrease(const Item& item, const Value& prio) { 1622 int id = index[item]; 1623 int kd = nodes[id].parent; 1624 nodes[id].prio = prio; 1625 while (kd >= 0 && less(id, kd)) { 1626 nodes[kd].prio = prio; 1627 nodes[kd].item = item; 1628 kd = nodes[kd].parent; 1629 } 1630 } 1631 1632 /// \brief Gives back the minimum priority of the class. 1633 /// 1634 /// \return Gives back the minimum priority of the class. 1635 const Value& classPrio(int cls) const { 1636 return nodes[~(classes[cls].parent)].prio; 1637 } 1638 1639 /// \brief Gives back the minimum priority item of the class. 1640 /// 1641 /// \return Gives back the minimum priority item of the class. 1642 const Item& classTop(int cls) const { 1643 return nodes[~(classes[cls].parent)].item; 1644 } 1645 1646 /// \brief Gives back a representant item of the class. 1647 /// 1648 /// The representant is indpendent from the priorities of the 1649 /// items. 1650 /// \return Gives back a representant item of the class. 1651 const Item& classRep(int id) const { 1652 int parent = classes[id].parent; 1653 return nodes[parent >= 0 ? classes[id].depth : leftNode(id)].item; 1654 } 1655 1656 /// \brief Lemon style iterator for the items of a class. 1657 /// 1658 /// ClassIt is a lemon style iterator for the components. It iterates 1659 /// on the items of a class. By example if you want to iterate on 1660 /// each items of each classes then you may write the next code. 1661 ///\code 1662 /// for (ClassIt cit(huf); cit != INVALID; ++cit) { 1663 /// std::cout << "Class: "; 1664 /// for (ItemIt iit(huf, cit); iit != INVALID; ++iit) { 1665 /// std::cout << toString(iit) << ' ' << std::endl; 1666 /// } 1667 /// std::cout << std::endl; 1668 /// } 1669 ///\endcode 1670 class ItemIt { 1671 private: 1672 1673 const HeapUnionFind* _huf; 1674 int _id, _lid; 1675 1676 public: 1677 1678 /// \brief Default constructor 1679 /// 1680 /// Default constructor 1681 ItemIt() {} 1682 1683 ItemIt(const HeapUnionFind& huf, int cls) : _huf(&huf) { 1684 int id = cls; 1685 int parent = _huf>classes[id].parent; 1686 if (parent >= 0) { 1687 _id = _huf>classes[id].depth; 1688 if (_huf>classes[id].next != 1) { 1689 _lid = _huf>classes[_huf>classes[id].next].depth; 1690 } else { 1691 _lid = 1; 1692 } 1693 } else { 1694 _id = _huf>leftNode(id); 1695 _lid = 1; 1696 } 1697 } 1698 1699 /// \brief Increment operator 1700 /// 1701 /// It steps to the next item in the class. 1702 ItemIt& operator++() { 1703 _id = _huf>nextNode(_id); 1704 return *this; 1705 } 1706 1707 /// \brief Conversion operator 1708 /// 1709 /// It converts the iterator to the current item. 1710 operator const Item&() const { 1711 return _huf>nodes[_id].item; 1712 } 1713 1714 /// \brief Equality operator 1715 /// 1716 /// Equality operator 1717 bool operator==(const ItemIt& i) { 1718 return i._id == _id; 1719 } 1720 1721 /// \brief Inequality operator 1722 /// 1723 /// Inequality operator 1724 bool operator!=(const ItemIt& i) { 1725 return i._id != _id; 1726 } 1727 1728 /// \brief Equality operator 1729 /// 1730 /// Equality operator 1731 bool operator==(Invalid) { 1732 return _id == _lid; 1733 } 1734 1735 /// \brief Inequality operator 1736 /// 1737 /// Inequality operator 1738 bool operator!=(Invalid) { 1739 return _id != _lid; 1740 } 1741 1742 }; 1743 1744 /// \brief Class iterator 1745 /// 1746 /// The iterator stores 1747 class ClassIt { 1748 private: 1749 1750 const HeapUnionFind* _huf; 1751 int _id; 1752 1753 public: 1754 1755 ClassIt(const HeapUnionFind& huf) 1756 : _huf(&huf), _id(huf.first_class) {} 1757 1758 ClassIt(const HeapUnionFind& huf, int cls) 1759 : _huf(&huf), _id(huf.classes[cls].left) {} 1760 1761 ClassIt(Invalid) : _huf(0), _id(1) {} 1762 1763 const ClassIt& operator++() { 1764 _id = _huf>classes[_id].next; 1765 return *this; 1766 } 1767 1768 /// \brief Equality operator 1769 /// 1770 /// Equality operator 1771 bool operator==(const ClassIt& i) { 1772 return i._id == _id; 1773 } 1774 1775 /// \brief Inequality operator 1776 /// 1777 /// Inequality operator 1778 bool operator!=(const ClassIt& i) { 1779 return i._id != _id; 1780 } 1781 1782 operator int() const { 1783 return _id; 1784 } 1785 1786 }; 1787 1788 }; 1789 927 1790 //! @} 928 1791
Note: See TracChangeset
for help on using the changeset viewer.