Changeset 2548:a3ba22ebccc6 in lemon-0.x for lemon
- Timestamp:
- 12/28/07 12:00:51 (16 years ago)
- Branch:
- default
- Phase:
- public
- Convert:
- svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3425
- Location:
- lemon
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
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 Union-Find 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 Union-Find 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 union-find 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 union-find. 1276 /// 1277 /// Constructs the union-find. 1278 /// \brief _index The index map of the union-find. 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 one-item-heap. 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 union-find 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 union-find 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.