Changeset 947:0513ccfea967 in lemon for lemon
 Timestamp:
 09/20/09 21:38:24 (12 years ago)
 Branch:
 default
 Phase:
 public
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/matching.h
r698 r947 17 17 */ 18 18 19 #ifndef LEMON_MA X_MATCHING_H20 #define LEMON_MA X_MATCHING_H19 #ifndef LEMON_MATCHING_H 20 #define LEMON_MATCHING_H 21 21 22 22 #include <vector> … … 42 42 /// This class implements Edmonds' alternating forest matching algorithm 43 43 /// for finding a maximum cardinality matching in a general undirected graph. 44 /// It can be started from an arbitrary initial matching 44 /// It can be started from an arbitrary initial matching 45 45 /// (the default is the empty one). 46 46 /// … … 70 70 ///\brief Status constants for GallaiEdmonds decomposition. 71 71 /// 72 ///These constants are used for indicating the GallaiEdmonds 72 ///These constants are used for indicating the GallaiEdmonds 73 73 ///decomposition of a graph. The nodes with status \c EVEN (or \c D) 74 74 ///induce a subgraph with factorcritical components, the nodes with 75 75 ///status \c ODD (or \c A) form the canonical barrier, and the nodes 76 ///with status \c MATCHED (or \c C) induce a subgraph having a 76 ///with status \c MATCHED (or \c C) induce a subgraph having a 77 77 ///perfect matching. 78 78 enum Status { … … 513 513 } 514 514 515 /// \brief Start Edmonds' algorithm with a heuristic improvement 515 /// \brief Start Edmonds' algorithm with a heuristic improvement 516 516 /// for dense graphs 517 517 /// … … 535 535 /// \brief Run Edmonds' algorithm 536 536 /// 537 /// This function runs Edmonds' algorithm. An additional heuristic of 538 /// postponing shrinks is used for relatively dense graphs 537 /// This function runs Edmonds' algorithm. An additional heuristic of 538 /// postponing shrinks is used for relatively dense graphs 539 539 /// (for which <tt>m>=2*n</tt> holds). 540 540 void run() { … … 557 557 /// \brief Return the size (cardinality) of the matching. 558 558 /// 559 /// This function returns the size (cardinality) of the current matching. 559 /// This function returns the size (cardinality) of the current matching. 560 560 /// After run() it returns the size of the maximum matching in the graph. 561 561 int matchingSize() const { … … 571 571 /// \brief Return \c true if the given edge is in the matching. 572 572 /// 573 /// This function returns \c true if the given edge is in the current 573 /// This function returns \c true if the given edge is in the current 574 574 /// matching. 575 575 bool matching(const Edge& edge) const { … … 580 580 /// 581 581 /// This function returns the matching arc (or edge) incident to the 582 /// given node in the current matching or \c INVALID if the node is 582 /// given node in the current matching or \c INVALID if the node is 583 583 /// not covered by the matching. 584 584 Arc matching(const Node& n) const { … … 596 596 /// \brief Return the mate of the given node. 597 597 /// 598 /// This function returns the mate of the given node in the current 598 /// This function returns the mate of the given node in the current 599 599 /// matching or \c INVALID if the node is not covered by the matching. 600 600 Node mate(const Node& n) const { … … 606 606 607 607 /// \name Dual Solution 608 /// Functions to get the dual solution, i.e. the GallaiEdmonds 608 /// Functions to get the dual solution, i.e. the GallaiEdmonds 609 609 /// decomposition. 610 610 … … 649 649 /// \f$O(nm\log n)\f$ time complexity. 650 650 /// 651 /// The maximum weighted matching problem is to find a subset of the 652 /// edges in an undirected graph with maximum overall weight for which 651 /// The maximum weighted matching problem is to find a subset of the 652 /// edges in an undirected graph with maximum overall weight for which 653 653 /// each node has at most one incident edge. 654 654 /// It can be formulated with the following linear program. … … 674 674 \frac{\vert B \vert  1}{2}z_B\f] */ 675 675 /// 676 /// The algorithm can be executed with the run() function. 676 /// The algorithm can be executed with the run() function. 677 677 /// After it the matching (the primal solution) and the dual solution 678 /// can be obtained using the query functions and the 679 /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, 680 /// which is able to iterate on the nodes of a blossom. 678 /// can be obtained using the query functions and the 679 /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, 680 /// which is able to iterate on the nodes of a blossom. 681 681 /// If the value type is integer, then the dual solution is multiplied 682 682 /// by \ref MaxWeightedMatching::dualScale "4". 683 683 /// 684 684 /// \tparam GR The undirected graph type the algorithm runs on. 685 /// \tparam WM The type edge weight map. The default type is 685 /// \tparam WM The type edge weight map. The default type is 686 686 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>". 687 687 #ifdef DOXYGEN … … 746 746 747 747 enum Status { 748 EVEN = 1, MATCHED = 0, ODD = 1 , UNMATCHED = 2748 EVEN = 1, MATCHED = 0, ODD = 1 749 749 }; 750 750 … … 845 845 846 846 void destroyStructures() { 847 _node_num = countNodes(_graph);848 _blossom_num = _node_num * 3 / 2;849 850 847 if (_matching) { 851 848 delete _matching; … … 922 919 if (_delta3>state(e) != _delta3>IN_HEAP && blossom != vb) { 923 920 _delta3>push(e, rw / 2); 924 }925 } else if ((*_blossom_data)[vb].status == UNMATCHED) {926 if (_delta3>state(e) != _delta3>IN_HEAP) {927 _delta3>push(e, rw);928 921 } 929 922 } else { … … 950 943 (*_blossom_data)[vb].offset); 951 944 } else if ((*_delta2)[vb] > _blossom_set>classPrio(vb)  952 (*_blossom_data)[vb].offset) {945 (*_blossom_data)[vb].offset) { 953 946 _delta2>decrease(vb, _blossom_set>classPrio(vb)  954 947 (*_blossom_data)[vb].offset); … … 969 962 if (!_blossom_set>trivial(blossom)) { 970 963 _delta4>push(blossom, (*_blossom_data)[blossom].pot / 2 + 971 (*_blossom_data)[blossom].offset);964 (*_blossom_data)[blossom].offset); 972 965 } 973 966 } … … 1037 1030 } 1038 1031 } 1039 1040 } else if ((*_blossom_data)[vb].status == UNMATCHED) {1041 if (_delta3>state(e) == _delta3>IN_HEAP) {1042 _delta3>erase(e);1043 }1044 1032 } else { 1045 1033 … … 1078 1066 std::numeric_limits<Value>::max()) { 1079 1067 _delta2>push(blossom, _blossom_set>classPrio(blossom)  1080 1068 (*_blossom_data)[blossom].offset); 1081 1069 } 1082 1070 … … 1117 1105 if (_delta3>state(e) != _delta3>IN_HEAP && blossom != vb) { 1118 1106 _delta3>push(e, rw / 2); 1119 }1120 } else if ((*_blossom_data)[vb].status == UNMATCHED) {1121 if (_delta3>state(e) != _delta3>IN_HEAP) {1122 _delta3>push(e, rw);1123 1107 } 1124 1108 } else { … … 1158 1142 } 1159 1143 1160 1161 void matchedToUnmatched(int blossom) {1162 if (_delta2>state(blossom) == _delta2>IN_HEAP) {1163 _delta2>erase(blossom);1164 }1165 1166 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);1167 n != INVALID; ++n) {1168 int ni = (*_node_index)[n];1169 1170 _blossom_set>increase(n, std::numeric_limits<Value>::max());1171 1172 (*_node_data)[ni].heap.clear();1173 (*_node_data)[ni].heap_index.clear();1174 1175 for (OutArcIt e(_graph, n); e != INVALID; ++e) {1176 Node v = _graph.target(e);1177 int vb = _blossom_set>find(v);1178 int vi = (*_node_index)[v];1179 1180 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot 1181 dualScale * _weight[e];1182 1183 if ((*_blossom_data)[vb].status == EVEN) {1184 if (_delta3>state(e) != _delta3>IN_HEAP) {1185 _delta3>push(e, rw);1186 }1187 }1188 }1189 }1190 }1191 1192 void unmatchedToMatched(int blossom) {1193 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);1194 n != INVALID; ++n) {1195 int ni = (*_node_index)[n];1196 1197 for (InArcIt e(_graph, n); e != INVALID; ++e) {1198 Node v = _graph.source(e);1199 int vb = _blossom_set>find(v);1200 int vi = (*_node_index)[v];1201 1202 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot 1203 dualScale * _weight[e];1204 1205 if (vb == blossom) {1206 if (_delta3>state(e) == _delta3>IN_HEAP) {1207 _delta3>erase(e);1208 }1209 } else if ((*_blossom_data)[vb].status == EVEN) {1210 1211 if (_delta3>state(e) == _delta3>IN_HEAP) {1212 _delta3>erase(e);1213 }1214 1215 int vt = _tree_set>find(vb);1216 1217 Arc r = _graph.oppositeArc(e);1218 1219 typename std::map<int, Arc>::iterator it =1220 (*_node_data)[ni].heap_index.find(vt);1221 1222 if (it != (*_node_data)[ni].heap_index.end()) {1223 if ((*_node_data)[ni].heap[it>second] > rw) {1224 (*_node_data)[ni].heap.replace(it>second, r);1225 (*_node_data)[ni].heap.decrease(r, rw);1226 it>second = r;1227 }1228 } else {1229 (*_node_data)[ni].heap.push(r, rw);1230 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));1231 }1232 1233 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {1234 _blossom_set>decrease(n, (*_node_data)[ni].heap.prio());1235 1236 if (_delta2>state(blossom) != _delta2>IN_HEAP) {1237 _delta2>push(blossom, _blossom_set>classPrio(blossom) 1238 (*_blossom_data)[blossom].offset);1239 } else if ((*_delta2)[blossom] > _blossom_set>classPrio(blossom)1240 (*_blossom_data)[blossom].offset){1241 _delta2>decrease(blossom, _blossom_set>classPrio(blossom) 1242 (*_blossom_data)[blossom].offset);1243 }1244 }1245 1246 } else if ((*_blossom_data)[vb].status == UNMATCHED) {1247 if (_delta3>state(e) == _delta3>IN_HEAP) {1248 _delta3>erase(e);1249 }1250 }1251 }1252 }1253 }1254 1255 1144 void alternatePath(int even, int tree) { 1256 1145 int odd; … … 1295 1184 destroyTree(tree); 1296 1185 1297 (*_blossom_data)[blossom].status = UNMATCHED;1298 1186 (*_blossom_data)[blossom].base = node; 1299 matchedToUnmatched(blossom); 1300 } 1301 1187 (*_blossom_data)[blossom].next = INVALID; 1188 } 1302 1189 1303 1190 void augmentOnEdge(const Edge& edge) { … … 1306 1193 int right = _blossom_set>find(_graph.v(edge)); 1307 1194 1308 if ((*_blossom_data)[left].status == EVEN) { 1309 int left_tree = _tree_set>find(left); 1310 alternatePath(left, left_tree); 1311 destroyTree(left_tree); 1312 } else { 1313 (*_blossom_data)[left].status = MATCHED; 1314 unmatchedToMatched(left); 1315 } 1316 1317 if ((*_blossom_data)[right].status == EVEN) { 1318 int right_tree = _tree_set>find(right); 1319 alternatePath(right, right_tree); 1320 destroyTree(right_tree); 1321 } else { 1322 (*_blossom_data)[right].status = MATCHED; 1323 unmatchedToMatched(right); 1324 } 1195 int left_tree = _tree_set>find(left); 1196 alternatePath(left, left_tree); 1197 destroyTree(left_tree); 1198 1199 int right_tree = _tree_set>find(right); 1200 alternatePath(right, right_tree); 1201 destroyTree(right_tree); 1325 1202 1326 1203 (*_blossom_data)[left].next = _graph.direct(edge, true); 1327 1204 (*_blossom_data)[right].next = _graph.direct(edge, false); 1205 } 1206 1207 void augmentOnArc(const Arc& arc) { 1208 1209 int left = _blossom_set>find(_graph.source(arc)); 1210 int right = _blossom_set>find(_graph.target(arc)); 1211 1212 (*_blossom_data)[left].status = MATCHED; 1213 1214 int right_tree = _tree_set>find(right); 1215 alternatePath(right, right_tree); 1216 destroyTree(right_tree); 1217 1218 (*_blossom_data)[left].next = arc; 1219 (*_blossom_data)[right].next = _graph.oppositeArc(arc); 1328 1220 } 1329 1221 … … 1530 1422 (*_blossom_data)[sb].pred = pred; 1531 1423 (*_blossom_data)[sb].next = 1532 1424 _graph.oppositeArc((*_blossom_data)[tb].next); 1533 1425 1534 1426 pred = (*_blossom_data)[ub].next; … … 1630 1522 1631 1523 for (int i = 0; i < int(blossoms.size()); ++i) { 1632 if ((*_blossom_data)[blossoms[i]]. status == MATCHED) {1524 if ((*_blossom_data)[blossoms[i]].next != INVALID) { 1633 1525 1634 1526 Value offset = (*_blossom_data)[blossoms[i]].offset; … … 1758 1650 _delta4>prio() : std::numeric_limits<Value>::max(); 1759 1651 1760 _delta_sum = d1; OpType ot = D1; 1652 _delta_sum = d3; OpType ot = D3; 1653 if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; } 1761 1654 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } 1762 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }1763 1655 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } 1764 1765 1656 1766 1657 switch (ot) { … … 1776 1667 int blossom = _delta2>top(); 1777 1668 Node n = _blossom_set>classTop(blossom); 1778 Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); 1779 extendOnArc(e); 1669 Arc a = (*_node_data)[(*_node_index)[n]].heap.top(); 1670 if ((*_blossom_data)[blossom].next == INVALID) { 1671 augmentOnArc(a); 1672 unmatched; 1673 } else { 1674 extendOnArc(a); 1675 } 1780 1676 } 1781 1677 break; … … 1790 1686 _delta3>pop(); 1791 1687 } else { 1792 int left_tree; 1793 if ((*_blossom_data)[left_blossom].status == EVEN) { 1794 left_tree = _tree_set>find(left_blossom); 1795 } else { 1796 left_tree = 1; 1797 ++unmatched; 1798 } 1799 int right_tree; 1800 if ((*_blossom_data)[right_blossom].status == EVEN) { 1801 right_tree = _tree_set>find(right_blossom); 1802 } else { 1803 right_tree = 1; 1804 ++unmatched; 1805 } 1688 int left_tree = _tree_set>find(left_blossom); 1689 int right_tree = _tree_set>find(right_blossom); 1806 1690 1807 1691 if (left_tree == right_tree) { … … 1838 1722 1839 1723 /// \name Primal Solution 1840 /// Functions to get the primal solution, i.e. the maximum weighted 1724 /// Functions to get the primal solution, i.e. the maximum weighted 1841 1725 /// matching.\n 1842 1726 /// Either \ref run() or \ref start() function should be called before … … 1857 1741 } 1858 1742 } 1859 return sum / =2;1743 return sum / 2; 1860 1744 } 1861 1745 … … 1877 1761 /// \brief Return \c true if the given edge is in the matching. 1878 1762 /// 1879 /// This function returns \c true if the given edge is in the found 1763 /// This function returns \c true if the given edge is in the found 1880 1764 /// matching. 1881 1765 /// … … 1888 1772 /// 1889 1773 /// This function returns the matching arc (or edge) incident to the 1890 /// given node in the found matching or \c INVALID if the node is 1774 /// given node in the found matching or \c INVALID if the node is 1891 1775 /// not covered by the matching. 1892 1776 /// … … 1906 1790 /// \brief Return the mate of the given node. 1907 1791 /// 1908 /// This function returns the mate of the given node in the found 1792 /// This function returns the mate of the given node in the found 1909 1793 /// matching or \c INVALID if the node is not covered by the matching. 1910 1794 /// … … 1926 1810 /// \brief Return the value of the dual solution. 1927 1811 /// 1928 /// This function returns the value of the dual solution. 1929 /// It should be equal to the primal value scaled by \ref dualScale 1812 /// This function returns the value of the dual solution. 1813 /// It should be equal to the primal value scaled by \ref dualScale 1930 1814 /// "dual scale". 1931 1815 /// … … 1982 1866 /// \brief Iterator for obtaining the nodes of a blossom. 1983 1867 /// 1984 /// This class provides an iterator for obtaining the nodes of the 1868 /// This class provides an iterator for obtaining the nodes of the 1985 1869 /// given blossom. It lists a subset of the nodes. 1986 /// Before using this iterator, you must allocate a 1870 /// Before using this iterator, you must allocate a 1987 1871 /// MaxWeightedMatching class and execute it. 1988 1872 class BlossomIt { … … 1993 1877 /// Constructor to get the nodes of the given variable. 1994 1878 /// 1995 /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 1996 /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 1879 /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 1880 /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 1997 1881 /// called before initializing this iterator. 1998 1882 BlossomIt(const MaxWeightedMatching& algorithm, int variable) … … 2047 1931 /// \f$O(nm\log n)\f$ time complexity. 2048 1932 /// 2049 /// The maximum weighted perfect matching problem is to find a subset of 2050 /// the edges in an undirected graph with maximum overall weight for which 1933 /// The maximum weighted perfect matching problem is to find a subset of 1934 /// the edges in an undirected graph with maximum overall weight for which 2051 1935 /// each node has exactly one incident edge. 2052 1936 /// It can be formulated with the following linear program. … … 2071 1955 \frac{\vert B \vert  1}{2}z_B\f] */ 2072 1956 /// 2073 /// The algorithm can be executed with the run() function. 1957 /// The algorithm can be executed with the run() function. 2074 1958 /// After it the matching (the primal solution) and the dual solution 2075 /// can be obtained using the query functions and the 2076 /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 2077 /// which is able to iterate on the nodes of a blossom. 1959 /// can be obtained using the query functions and the 1960 /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 1961 /// which is able to iterate on the nodes of a blossom. 2078 1962 /// If the value type is integer, then the dual solution is multiplied 2079 1963 /// by \ref MaxWeightedMatching::dualScale "4". 2080 1964 /// 2081 1965 /// \tparam GR The undirected graph type the algorithm runs on. 2082 /// \tparam WM The type edge weight map. The default type is 1966 /// \tparam WM The type edge weight map. The default type is 2083 1967 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>". 2084 1968 #ifdef DOXYGEN … … 2234 2118 2235 2119 void destroyStructures() { 2236 _node_num = countNodes(_graph);2237 _blossom_num = _node_num * 3 / 2;2238 2239 2120 if (_matching) { 2240 2121 delete _matching; … … 2992 2873 _delta4>prio() : std::numeric_limits<Value>::max(); 2993 2874 2994 _delta_sum = d 2; OpType ot = D2;2995 if (d 3 < _delta_sum) { _delta_sum = d3; ot = D3; }2875 _delta_sum = d3; OpType ot = D3; 2876 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } 2996 2877 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } 2997 2878 … … 3056 2937 3057 2938 /// \name Primal Solution 3058 /// Functions to get the primal solution, i.e. the maximum weighted 2939 /// Functions to get the primal solution, i.e. the maximum weighted 3059 2940 /// perfect matching.\n 3060 2941 /// Either \ref run() or \ref start() function should be called before … … 3075 2956 } 3076 2957 } 3077 return sum / =2;2958 return sum / 2; 3078 2959 } 3079 2960 3080 2961 /// \brief Return \c true if the given edge is in the matching. 3081 2962 /// 3082 /// This function returns \c true if the given edge is in the found 2963 /// This function returns \c true if the given edge is in the found 3083 2964 /// matching. 3084 2965 /// … … 3091 2972 /// 3092 2973 /// This function returns the matching arc (or edge) incident to the 3093 /// given node in the found matching or \c INVALID if the node is 2974 /// given node in the found matching or \c INVALID if the node is 3094 2975 /// not covered by the matching. 3095 2976 /// … … 3109 2990 /// \brief Return the mate of the given node. 3110 2991 /// 3111 /// This function returns the mate of the given node in the found 2992 /// This function returns the mate of the given node in the found 3112 2993 /// matching or \c INVALID if the node is not covered by the matching. 3113 2994 /// … … 3128 3009 /// \brief Return the value of the dual solution. 3129 3010 /// 3130 /// This function returns the value of the dual solution. 3131 /// It should be equal to the primal value scaled by \ref dualScale 3011 /// This function returns the value of the dual solution. 3012 /// It should be equal to the primal value scaled by \ref dualScale 3132 3013 /// "dual scale". 3133 3014 /// … … 3184 3065 /// \brief Iterator for obtaining the nodes of a blossom. 3185 3066 /// 3186 /// This class provides an iterator for obtaining the nodes of the 3067 /// This class provides an iterator for obtaining the nodes of the 3187 3068 /// given blossom. It lists a subset of the nodes. 3188 /// Before using this iterator, you must allocate a 3069 /// Before using this iterator, you must allocate a 3189 3070 /// MaxWeightedPerfectMatching class and execute it. 3190 3071 class BlossomIt { … … 3195 3076 /// Constructor to get the nodes of the given variable. 3196 3077 /// 3197 /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 3198 /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 3078 /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 3079 /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 3199 3080 /// must be called before initializing this iterator. 3200 3081 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) … … 3242 3123 } //END OF NAMESPACE LEMON 3243 3124 3244 #endif //LEMON_MA X_MATCHING_H3125 #endif //LEMON_MATCHING_H
Note: See TracChangeset
for help on using the changeset viewer.