Changeset 954:07ec2b52e53d in lemon for lemon/matching.h
 Timestamp:
 03/17/10 10:29:57 (10 years ago)
 Branch:
 default
 Parents:
 953:d8ea85825e02 (diff), 945:5b926cc36a4b (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.  Phase:
 public
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

lemon/matching.h
r945 r954 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> … … 29 29 #include <lemon/bin_heap.h> 30 30 #include <lemon/maps.h> 31 #include <lemon/fractional_matching.h> 31 32 32 33 ///\ingroup matching … … 42 43 /// This class implements Edmonds' alternating forest matching algorithm 43 44 /// for finding a maximum cardinality matching in a general undirected graph. 44 /// It can be started from an arbitrary initial matching 45 /// It can be started from an arbitrary initial matching 45 46 /// (the default is the empty one). 46 47 /// … … 70 71 ///\brief Status constants for GallaiEdmonds decomposition. 71 72 /// 72 ///These constants are used for indicating the GallaiEdmonds 73 ///These constants are used for indicating the GallaiEdmonds 73 74 ///decomposition of a graph. The nodes with status \c EVEN (or \c D) 74 75 ///induce a subgraph with factorcritical components, the nodes with 75 76 ///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 77 ///with status \c MATCHED (or \c C) induce a subgraph having a 77 78 ///perfect matching. 78 79 enum Status { … … 513 514 } 514 515 515 /// \brief Start Edmonds' algorithm with a heuristic improvement 516 /// \brief Start Edmonds' algorithm with a heuristic improvement 516 517 /// for dense graphs 517 518 /// … … 535 536 /// \brief Run Edmonds' algorithm 536 537 /// 537 /// This function runs Edmonds' algorithm. An additional heuristic of 538 /// postponing shrinks is used for relatively dense graphs 538 /// This function runs Edmonds' algorithm. An additional heuristic of 539 /// postponing shrinks is used for relatively dense graphs 539 540 /// (for which <tt>m>=2*n</tt> holds). 540 541 void run() { … … 557 558 /// \brief Return the size (cardinality) of the matching. 558 559 /// 559 /// This function returns the size (cardinality) of the current matching. 560 /// This function returns the size (cardinality) of the current matching. 560 561 /// After run() it returns the size of the maximum matching in the graph. 561 562 int matchingSize() const { … … 571 572 /// \brief Return \c true if the given edge is in the matching. 572 573 /// 573 /// This function returns \c true if the given edge is in the current 574 /// This function returns \c true if the given edge is in the current 574 575 /// matching. 575 576 bool matching(const Edge& edge) const { … … 580 581 /// 581 582 /// 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 583 /// given node in the current matching or \c INVALID if the node is 583 584 /// not covered by the matching. 584 585 Arc matching(const Node& n) const { … … 596 597 /// \brief Return the mate of the given node. 597 598 /// 598 /// This function returns the mate of the given node in the current 599 /// This function returns the mate of the given node in the current 599 600 /// matching or \c INVALID if the node is not covered by the matching. 600 601 Node mate(const Node& n) const { … … 606 607 607 608 /// \name Dual Solution 608 /// Functions to get the dual solution, i.e. the GallaiEdmonds 609 /// Functions to get the dual solution, i.e. the GallaiEdmonds 609 610 /// decomposition. 610 611 … … 649 650 /// \f$O(nm\log n)\f$ time complexity. 650 651 /// 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 652 /// The maximum weighted matching problem is to find a subset of the 653 /// edges in an undirected graph with maximum overall weight for which 653 654 /// each node has at most one incident edge. 654 655 /// It can be formulated with the following linear program. … … 674 675 \frac{\vert B \vert  1}{2}z_B\f] */ 675 676 /// 676 /// The algorithm can be executed with the run() function. 677 /// The algorithm can be executed with the run() function. 677 678 /// 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. 679 /// can be obtained using the query functions and the 680 /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, 681 /// which is able to iterate on the nodes of a blossom. 681 682 /// If the value type is integer, then the dual solution is multiplied 682 683 /// by \ref MaxWeightedMatching::dualScale "4". 683 684 /// 684 685 /// \tparam GR The undirected graph type the algorithm runs on. 685 /// \tparam WM The type edge weight map. The default type is 686 /// \tparam WM The type edge weight map. The default type is 686 687 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>". 687 688 #ifdef DOXYGEN … … 746 747 747 748 enum Status { 748 EVEN = 1, MATCHED = 0, ODD = 1 , UNMATCHED = 2749 EVEN = 1, MATCHED = 0, ODD = 1 749 750 }; 750 751 … … 798 799 799 800 Value _delta_sum; 801 int _unmatched; 802 803 typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching; 804 FractionalMatching *_fractional; 800 805 801 806 void createStructures() { … … 864 869 865 870 void destroyStructures() { 866 _node_num = countNodes(_graph);867 _blossom_num = _node_num * 3 / 2;868 869 871 if (_matching) { 870 872 delete _matching; … … 941 943 if (_delta3>state(e) != _delta3>IN_HEAP && blossom != vb) { 942 944 _delta3>push(e, rw / 2); 943 }944 } else if ((*_blossom_data)[vb].status == UNMATCHED) {945 if (_delta3>state(e) != _delta3>IN_HEAP) {946 _delta3>push(e, rw);947 945 } 948 946 } else { … … 969 967 (*_blossom_data)[vb].offset); 970 968 } else if ((*_delta2)[vb] > _blossom_set>classPrio(vb)  971 (*_blossom_data)[vb].offset) {969 (*_blossom_data)[vb].offset) { 972 970 _delta2>decrease(vb, _blossom_set>classPrio(vb)  973 971 (*_blossom_data)[vb].offset); … … 988 986 if (!_blossom_set>trivial(blossom)) { 989 987 _delta4>push(blossom, (*_blossom_data)[blossom].pot / 2 + 990 (*_blossom_data)[blossom].offset);988 (*_blossom_data)[blossom].offset); 991 989 } 992 990 } … … 1055 1053 } 1056 1054 } 1057 }1058 1059 } else if ((*_blossom_data)[vb].status == UNMATCHED) {1060 if (_delta3>state(e) == _delta3>IN_HEAP) {1061 _delta3>erase(e);1062 1055 } 1063 1056 } else { … … 1097 1090 std::numeric_limits<Value>::max()) { 1098 1091 _delta2>push(blossom, _blossom_set>classPrio(blossom)  1099 1092 (*_blossom_data)[blossom].offset); 1100 1093 } 1101 1094 … … 1136 1129 if (_delta3>state(e) != _delta3>IN_HEAP && blossom != vb) { 1137 1130 _delta3>push(e, rw / 2); 1138 }1139 } else if ((*_blossom_data)[vb].status == UNMATCHED) {1140 if (_delta3>state(e) != _delta3>IN_HEAP) {1141 _delta3>push(e, rw);1142 1131 } 1143 1132 } else { … … 1177 1166 } 1178 1167 1179 1180 void matchedToUnmatched(int blossom) {1181 if (_delta2>state(blossom) == _delta2>IN_HEAP) {1182 _delta2>erase(blossom);1183 }1184 1185 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);1186 n != INVALID; ++n) {1187 int ni = (*_node_index)[n];1188 1189 _blossom_set>increase(n, std::numeric_limits<Value>::max());1190 1191 (*_node_data)[ni].heap.clear();1192 (*_node_data)[ni].heap_index.clear();1193 1194 for (OutArcIt e(_graph, n); e != INVALID; ++e) {1195 Node v = _graph.target(e);1196 int vb = _blossom_set>find(v);1197 int vi = (*_node_index)[v];1198 1199 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot 1200 dualScale * _weight[e];1201 1202 if ((*_blossom_data)[vb].status == EVEN) {1203 if (_delta3>state(e) != _delta3>IN_HEAP) {1204 _delta3>push(e, rw);1205 }1206 }1207 }1208 }1209 }1210 1211 void unmatchedToMatched(int blossom) {1212 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);1213 n != INVALID; ++n) {1214 int ni = (*_node_index)[n];1215 1216 for (InArcIt e(_graph, n); e != INVALID; ++e) {1217 Node v = _graph.source(e);1218 int vb = _blossom_set>find(v);1219 int vi = (*_node_index)[v];1220 1221 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot 1222 dualScale * _weight[e];1223 1224 if (vb == blossom) {1225 if (_delta3>state(e) == _delta3>IN_HEAP) {1226 _delta3>erase(e);1227 }1228 } else if ((*_blossom_data)[vb].status == EVEN) {1229 1230 if (_delta3>state(e) == _delta3>IN_HEAP) {1231 _delta3>erase(e);1232 }1233 1234 int vt = _tree_set>find(vb);1235 1236 Arc r = _graph.oppositeArc(e);1237 1238 typename std::map<int, Arc>::iterator it =1239 (*_node_data)[ni].heap_index.find(vt);1240 1241 if (it != (*_node_data)[ni].heap_index.end()) {1242 if ((*_node_data)[ni].heap[it>second] > rw) {1243 (*_node_data)[ni].heap.replace(it>second, r);1244 (*_node_data)[ni].heap.decrease(r, rw);1245 it>second = r;1246 }1247 } else {1248 (*_node_data)[ni].heap.push(r, rw);1249 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));1250 }1251 1252 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {1253 _blossom_set>decrease(n, (*_node_data)[ni].heap.prio());1254 1255 if (_delta2>state(blossom) != _delta2>IN_HEAP) {1256 _delta2>push(blossom, _blossom_set>classPrio(blossom) 1257 (*_blossom_data)[blossom].offset);1258 } else if ((*_delta2)[blossom] > _blossom_set>classPrio(blossom)1259 (*_blossom_data)[blossom].offset){1260 _delta2>decrease(blossom, _blossom_set>classPrio(blossom) 1261 (*_blossom_data)[blossom].offset);1262 }1263 }1264 1265 } else if ((*_blossom_data)[vb].status == UNMATCHED) {1266 if (_delta3>state(e) == _delta3>IN_HEAP) {1267 _delta3>erase(e);1268 }1269 }1270 }1271 }1272 }1273 1274 1168 void alternatePath(int even, int tree) { 1275 1169 int odd; … … 1314 1208 destroyTree(tree); 1315 1209 1316 (*_blossom_data)[blossom].status = UNMATCHED;1317 1210 (*_blossom_data)[blossom].base = node; 1318 matchedToUnmatched(blossom); 1319 } 1320 1211 (*_blossom_data)[blossom].next = INVALID; 1212 } 1321 1213 1322 1214 void augmentOnEdge(const Edge& edge) { … … 1325 1217 int right = _blossom_set>find(_graph.v(edge)); 1326 1218 1327 if ((*_blossom_data)[left].status == EVEN) { 1328 int left_tree = _tree_set>find(left); 1329 alternatePath(left, left_tree); 1330 destroyTree(left_tree); 1331 } else { 1332 (*_blossom_data)[left].status = MATCHED; 1333 unmatchedToMatched(left); 1334 } 1335 1336 if ((*_blossom_data)[right].status == EVEN) { 1337 int right_tree = _tree_set>find(right); 1338 alternatePath(right, right_tree); 1339 destroyTree(right_tree); 1340 } else { 1341 (*_blossom_data)[right].status = MATCHED; 1342 unmatchedToMatched(right); 1343 } 1219 int left_tree = _tree_set>find(left); 1220 alternatePath(left, left_tree); 1221 destroyTree(left_tree); 1222 1223 int right_tree = _tree_set>find(right); 1224 alternatePath(right, right_tree); 1225 destroyTree(right_tree); 1344 1226 1345 1227 (*_blossom_data)[left].next = _graph.direct(edge, true); 1346 1228 (*_blossom_data)[right].next = _graph.direct(edge, false); 1229 } 1230 1231 void augmentOnArc(const Arc& arc) { 1232 1233 int left = _blossom_set>find(_graph.source(arc)); 1234 int right = _blossom_set>find(_graph.target(arc)); 1235 1236 (*_blossom_data)[left].status = MATCHED; 1237 1238 int right_tree = _tree_set>find(right); 1239 alternatePath(right, right_tree); 1240 destroyTree(right_tree); 1241 1242 (*_blossom_data)[left].next = arc; 1243 (*_blossom_data)[right].next = _graph.oppositeArc(arc); 1347 1244 } 1348 1245 … … 1549 1446 (*_blossom_data)[sb].pred = pred; 1550 1447 (*_blossom_data)[sb].next = 1551 1448 _graph.oppositeArc((*_blossom_data)[tb].next); 1552 1449 1553 1450 pred = (*_blossom_data)[ub].next; … … 1649 1546 1650 1547 for (int i = 0; i < int(blossoms.size()); ++i) { 1651 if ((*_blossom_data)[blossoms[i]]. status == MATCHED) {1548 if ((*_blossom_data)[blossoms[i]].next != INVALID) { 1652 1549 1653 1550 Value offset = (*_blossom_data)[blossoms[i]].offset; … … 1687 1584 _delta4_index(0), _delta4(0), 1688 1585 1689 _delta_sum() {} 1586 _delta_sum(), _unmatched(0), 1587 1588 _fractional(0) 1589 {} 1690 1590 1691 1591 ~MaxWeightedMatching() { 1692 1592 destroyStructures(); 1593 if (_fractional) { 1594 delete _fractional; 1595 } 1693 1596 } 1694 1597 … … 1722 1625 } 1723 1626 1627 _unmatched = _node_num; 1628 1724 1629 _delta1>clear(); 1725 1630 _delta2>clear(); … … 1765 1670 } 1766 1671 1672 /// \brief Initialize the algorithm with fractional matching 1673 /// 1674 /// This function initializes the algorithm with a fractional 1675 /// matching. This initialization is also called jumpstart heuristic. 1676 void fractionalInit() { 1677 createStructures(); 1678 1679 if (_fractional == 0) { 1680 _fractional = new FractionalMatching(_graph, _weight, false); 1681 } 1682 _fractional>run(); 1683 1684 for (ArcIt e(_graph); e != INVALID; ++e) { 1685 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; 1686 } 1687 for (NodeIt n(_graph); n != INVALID; ++n) { 1688 (*_delta1_index)[n] = _delta1>PRE_HEAP; 1689 } 1690 for (EdgeIt e(_graph); e != INVALID; ++e) { 1691 (*_delta3_index)[e] = _delta3>PRE_HEAP; 1692 } 1693 for (int i = 0; i < _blossom_num; ++i) { 1694 (*_delta2_index)[i] = _delta2>PRE_HEAP; 1695 (*_delta4_index)[i] = _delta4>PRE_HEAP; 1696 } 1697 1698 _unmatched = 0; 1699 1700 int index = 0; 1701 for (NodeIt n(_graph); n != INVALID; ++n) { 1702 Value pot = _fractional>nodeValue(n); 1703 (*_node_index)[n] = index; 1704 (*_node_data)[index].pot = pot; 1705 int blossom = 1706 _blossom_set>insert(n, std::numeric_limits<Value>::max()); 1707 1708 (*_blossom_data)[blossom].status = MATCHED; 1709 (*_blossom_data)[blossom].pred = INVALID; 1710 (*_blossom_data)[blossom].next = _fractional>matching(n); 1711 if (_fractional>matching(n) == INVALID) { 1712 (*_blossom_data)[blossom].base = n; 1713 } 1714 (*_blossom_data)[blossom].pot = 0; 1715 (*_blossom_data)[blossom].offset = 0; 1716 ++index; 1717 } 1718 1719 typename Graph::template NodeMap<bool> processed(_graph, false); 1720 for (NodeIt n(_graph); n != INVALID; ++n) { 1721 if (processed[n]) continue; 1722 processed[n] = true; 1723 if (_fractional>matching(n) == INVALID) continue; 1724 int num = 1; 1725 Node v = _graph.target(_fractional>matching(n)); 1726 while (n != v) { 1727 processed[v] = true; 1728 v = _graph.target(_fractional>matching(v)); 1729 ++num; 1730 } 1731 1732 if (num % 2 == 1) { 1733 std::vector<int> subblossoms(num); 1734 1735 subblossoms[num] = _blossom_set>find(n); 1736 _delta1>push(n, _fractional>nodeValue(n)); 1737 v = _graph.target(_fractional>matching(n)); 1738 while (n != v) { 1739 subblossoms[num] = _blossom_set>find(v); 1740 _delta1>push(v, _fractional>nodeValue(v)); 1741 v = _graph.target(_fractional>matching(v)); 1742 } 1743 1744 int surface = 1745 _blossom_set>join(subblossoms.begin(), subblossoms.end()); 1746 (*_blossom_data)[surface].status = EVEN; 1747 (*_blossom_data)[surface].pred = INVALID; 1748 (*_blossom_data)[surface].next = INVALID; 1749 (*_blossom_data)[surface].pot = 0; 1750 (*_blossom_data)[surface].offset = 0; 1751 1752 _tree_set>insert(surface); 1753 ++_unmatched; 1754 } 1755 } 1756 1757 for (EdgeIt e(_graph); e != INVALID; ++e) { 1758 int si = (*_node_index)[_graph.u(e)]; 1759 int sb = _blossom_set>find(_graph.u(e)); 1760 int ti = (*_node_index)[_graph.v(e)]; 1761 int tb = _blossom_set>find(_graph.v(e)); 1762 if ((*_blossom_data)[sb].status == EVEN && 1763 (*_blossom_data)[tb].status == EVEN && sb != tb) { 1764 _delta3>push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot  1765 dualScale * _weight[e]) / 2); 1766 } 1767 } 1768 1769 for (NodeIt n(_graph); n != INVALID; ++n) { 1770 int nb = _blossom_set>find(n); 1771 if ((*_blossom_data)[nb].status != MATCHED) continue; 1772 int ni = (*_node_index)[n]; 1773 1774 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 1775 Node v = _graph.target(e); 1776 int vb = _blossom_set>find(v); 1777 int vi = (*_node_index)[v]; 1778 1779 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot  1780 dualScale * _weight[e]; 1781 1782 if ((*_blossom_data)[vb].status == EVEN) { 1783 1784 int vt = _tree_set>find(vb); 1785 1786 typename std::map<int, Arc>::iterator it = 1787 (*_node_data)[ni].heap_index.find(vt); 1788 1789 if (it != (*_node_data)[ni].heap_index.end()) { 1790 if ((*_node_data)[ni].heap[it>second] > rw) { 1791 (*_node_data)[ni].heap.replace(it>second, e); 1792 (*_node_data)[ni].heap.decrease(e, rw); 1793 it>second = e; 1794 } 1795 } else { 1796 (*_node_data)[ni].heap.push(e, rw); 1797 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); 1798 } 1799 } 1800 } 1801 1802 if (!(*_node_data)[ni].heap.empty()) { 1803 _blossom_set>decrease(n, (*_node_data)[ni].heap.prio()); 1804 _delta2>push(nb, _blossom_set>classPrio(nb)); 1805 } 1806 } 1807 } 1808 1767 1809 /// \brief Start the algorithm 1768 1810 /// 1769 1811 /// This function starts the algorithm. 1770 1812 /// 1771 /// \pre \ref init() must be called before using this function. 1813 /// \pre \ref init() or \ref fractionalInit() must be called 1814 /// before using this function. 1772 1815 void start() { 1773 1816 enum OpType { … … 1775 1818 }; 1776 1819 1777 int unmatched = _node_num; 1778 while (unmatched > 0) { 1820 while (_unmatched > 0) { 1779 1821 Value d1 = !_delta1>empty() ? 1780 1822 _delta1>prio() : std::numeric_limits<Value>::max(); … … 1789 1831 _delta4>prio() : std::numeric_limits<Value>::max(); 1790 1832 1791 _delta_sum = d1; OpType ot = D1; 1833 _delta_sum = d3; OpType ot = D3; 1834 if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; } 1792 1835 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } 1793 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }1794 1836 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } 1795 1796 1837 1797 1838 switch (ot) { … … 1800 1841 Node n = _delta1>top(); 1801 1842 unmatchNode(n); 1802  unmatched;1843 _unmatched; 1803 1844 } 1804 1845 break; … … 1807 1848 int blossom = _delta2>top(); 1808 1849 Node n = _blossom_set>classTop(blossom); 1809 Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); 1810 extendOnArc(e); 1850 Arc a = (*_node_data)[(*_node_index)[n]].heap.top(); 1851 if ((*_blossom_data)[blossom].next == INVALID) { 1852 augmentOnArc(a); 1853 _unmatched; 1854 } else { 1855 extendOnArc(a); 1856 } 1811 1857 } 1812 1858 break; … … 1821 1867 _delta3>pop(); 1822 1868 } else { 1823 int left_tree; 1824 if ((*_blossom_data)[left_blossom].status == EVEN) { 1825 left_tree = _tree_set>find(left_blossom); 1826 } else { 1827 left_tree = 1; 1828 ++unmatched; 1829 } 1830 int right_tree; 1831 if ((*_blossom_data)[right_blossom].status == EVEN) { 1832 right_tree = _tree_set>find(right_blossom); 1833 } else { 1834 right_tree = 1; 1835 ++unmatched; 1836 } 1869 int left_tree = _tree_set>find(left_blossom); 1870 int right_tree = _tree_set>find(right_blossom); 1837 1871 1838 1872 if (left_tree == right_tree) { … … 1840 1874 } else { 1841 1875 augmentOnEdge(e); 1842 unmatched = 2;1876 _unmatched = 2; 1843 1877 } 1844 1878 } … … 1858 1892 /// \note mwm.run() is just a shortcut of the following code. 1859 1893 /// \code 1860 /// mwm. init();1894 /// mwm.fractionalInit(); 1861 1895 /// mwm.start(); 1862 1896 /// \endcode 1863 1897 void run() { 1864 init();1898 fractionalInit(); 1865 1899 start(); 1866 1900 } … … 1869 1903 1870 1904 /// \name Primal Solution 1871 /// Functions to get the primal solution, i.e. the maximum weighted 1905 /// Functions to get the primal solution, i.e. the maximum weighted 1872 1906 /// matching.\n 1873 1907 /// Either \ref run() or \ref start() function should be called before … … 1888 1922 } 1889 1923 } 1890 return sum / =2;1924 return sum / 2; 1891 1925 } 1892 1926 … … 1908 1942 /// \brief Return \c true if the given edge is in the matching. 1909 1943 /// 1910 /// This function returns \c true if the given edge is in the found 1944 /// This function returns \c true if the given edge is in the found 1911 1945 /// matching. 1912 1946 /// … … 1919 1953 /// 1920 1954 /// This function returns the matching arc (or edge) incident to the 1921 /// given node in the found matching or \c INVALID if the node is 1955 /// given node in the found matching or \c INVALID if the node is 1922 1956 /// not covered by the matching. 1923 1957 /// … … 1937 1971 /// \brief Return the mate of the given node. 1938 1972 /// 1939 /// This function returns the mate of the given node in the found 1973 /// This function returns the mate of the given node in the found 1940 1974 /// matching or \c INVALID if the node is not covered by the matching. 1941 1975 /// … … 1957 1991 /// \brief Return the value of the dual solution. 1958 1992 /// 1959 /// This function returns the value of the dual solution. 1960 /// It should be equal to the primal value scaled by \ref dualScale 1993 /// This function returns the value of the dual solution. 1994 /// It should be equal to the primal value scaled by \ref dualScale 1961 1995 /// "dual scale". 1962 1996 /// … … 2013 2047 /// \brief Iterator for obtaining the nodes of a blossom. 2014 2048 /// 2015 /// This class provides an iterator for obtaining the nodes of the 2049 /// This class provides an iterator for obtaining the nodes of the 2016 2050 /// given blossom. It lists a subset of the nodes. 2017 /// Before using this iterator, you must allocate a 2051 /// Before using this iterator, you must allocate a 2018 2052 /// MaxWeightedMatching class and execute it. 2019 2053 class BlossomIt { … … 2024 2058 /// Constructor to get the nodes of the given variable. 2025 2059 /// 2026 /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 2027 /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 2060 /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 2061 /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 2028 2062 /// called before initializing this iterator. 2029 2063 BlossomIt(const MaxWeightedMatching& algorithm, int variable) … … 2078 2112 /// \f$O(nm\log n)\f$ time complexity. 2079 2113 /// 2080 /// The maximum weighted perfect matching problem is to find a subset of 2081 /// the edges in an undirected graph with maximum overall weight for which 2114 /// The maximum weighted perfect matching problem is to find a subset of 2115 /// the edges in an undirected graph with maximum overall weight for which 2082 2116 /// each node has exactly one incident edge. 2083 2117 /// It can be formulated with the following linear program. … … 2102 2136 \frac{\vert B \vert  1}{2}z_B\f] */ 2103 2137 /// 2104 /// The algorithm can be executed with the run() function. 2138 /// The algorithm can be executed with the run() function. 2105 2139 /// After it the matching (the primal solution) and the dual solution 2106 /// can be obtained using the query functions and the 2107 /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 2108 /// which is able to iterate on the nodes of a blossom. 2140 /// can be obtained using the query functions and the 2141 /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 2142 /// which is able to iterate on the nodes of a blossom. 2109 2143 /// If the value type is integer, then the dual solution is multiplied 2110 2144 /// by \ref MaxWeightedMatching::dualScale "4". 2111 2145 /// 2112 2146 /// \tparam GR The undirected graph type the algorithm runs on. 2113 /// \tparam WM The type edge weight map. The default type is 2147 /// \tparam WM The type edge weight map. The default type is 2114 2148 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>". 2115 2149 #ifdef DOXYGEN … … 2222 2256 2223 2257 Value _delta_sum; 2258 int _unmatched; 2259 2260 typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap> 2261 FractionalMatching; 2262 FractionalMatching *_fractional; 2224 2263 2225 2264 void createStructures() { … … 2283 2322 2284 2323 void destroyStructures() { 2285 _node_num = countNodes(_graph);2286 _blossom_num = _node_num * 3 / 2;2287 2288 2324 if (_matching) { 2289 2325 delete _matching; … … 2958 2994 _delta4_index(0), _delta4(0), 2959 2995 2960 _delta_sum() {} 2996 _delta_sum(), _unmatched(0), 2997 2998 _fractional(0) 2999 {} 2961 3000 2962 3001 ~MaxWeightedPerfectMatching() { 2963 3002 destroyStructures(); 3003 if (_fractional) { 3004 delete _fractional; 3005 } 2964 3006 } 2965 3007 … … 2989 3031 (*_delta4_index)[i] = _delta4>PRE_HEAP; 2990 3032 } 3033 3034 _unmatched = _node_num; 2991 3035 2992 3036 _delta2>clear(); … … 3031 3075 } 3032 3076 3077 /// \brief Initialize the algorithm with fractional matching 3078 /// 3079 /// This function initializes the algorithm with a fractional 3080 /// matching. This initialization is also called jumpstart heuristic. 3081 void fractionalInit() { 3082 createStructures(); 3083 3084 if (_fractional == 0) { 3085 _fractional = new FractionalMatching(_graph, _weight, false); 3086 } 3087 if (!_fractional>run()) { 3088 _unmatched = 1; 3089 return; 3090 } 3091 3092 for (ArcIt e(_graph); e != INVALID; ++e) { 3093 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; 3094 } 3095 for (EdgeIt e(_graph); e != INVALID; ++e) { 3096 (*_delta3_index)[e] = _delta3>PRE_HEAP; 3097 } 3098 for (int i = 0; i < _blossom_num; ++i) { 3099 (*_delta2_index)[i] = _delta2>PRE_HEAP; 3100 (*_delta4_index)[i] = _delta4>PRE_HEAP; 3101 } 3102 3103 _unmatched = 0; 3104 3105 int index = 0; 3106 for (NodeIt n(_graph); n != INVALID; ++n) { 3107 Value pot = _fractional>nodeValue(n); 3108 (*_node_index)[n] = index; 3109 (*_node_data)[index].pot = pot; 3110 int blossom = 3111 _blossom_set>insert(n, std::numeric_limits<Value>::max()); 3112 3113 (*_blossom_data)[blossom].status = MATCHED; 3114 (*_blossom_data)[blossom].pred = INVALID; 3115 (*_blossom_data)[blossom].next = _fractional>matching(n); 3116 (*_blossom_data)[blossom].pot = 0; 3117 (*_blossom_data)[blossom].offset = 0; 3118 ++index; 3119 } 3120 3121 typename Graph::template NodeMap<bool> processed(_graph, false); 3122 for (NodeIt n(_graph); n != INVALID; ++n) { 3123 if (processed[n]) continue; 3124 processed[n] = true; 3125 if (_fractional>matching(n) == INVALID) continue; 3126 int num = 1; 3127 Node v = _graph.target(_fractional>matching(n)); 3128 while (n != v) { 3129 processed[v] = true; 3130 v = _graph.target(_fractional>matching(v)); 3131 ++num; 3132 } 3133 3134 if (num % 2 == 1) { 3135 std::vector<int> subblossoms(num); 3136 3137 subblossoms[num] = _blossom_set>find(n); 3138 v = _graph.target(_fractional>matching(n)); 3139 while (n != v) { 3140 subblossoms[num] = _blossom_set>find(v); 3141 v = _graph.target(_fractional>matching(v)); 3142 } 3143 3144 int surface = 3145 _blossom_set>join(subblossoms.begin(), subblossoms.end()); 3146 (*_blossom_data)[surface].status = EVEN; 3147 (*_blossom_data)[surface].pred = INVALID; 3148 (*_blossom_data)[surface].next = INVALID; 3149 (*_blossom_data)[surface].pot = 0; 3150 (*_blossom_data)[surface].offset = 0; 3151 3152 _tree_set>insert(surface); 3153 ++_unmatched; 3154 } 3155 } 3156 3157 for (EdgeIt e(_graph); e != INVALID; ++e) { 3158 int si = (*_node_index)[_graph.u(e)]; 3159 int sb = _blossom_set>find(_graph.u(e)); 3160 int ti = (*_node_index)[_graph.v(e)]; 3161 int tb = _blossom_set>find(_graph.v(e)); 3162 if ((*_blossom_data)[sb].status == EVEN && 3163 (*_blossom_data)[tb].status == EVEN && sb != tb) { 3164 _delta3>push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot  3165 dualScale * _weight[e]) / 2); 3166 } 3167 } 3168 3169 for (NodeIt n(_graph); n != INVALID; ++n) { 3170 int nb = _blossom_set>find(n); 3171 if ((*_blossom_data)[nb].status != MATCHED) continue; 3172 int ni = (*_node_index)[n]; 3173 3174 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 3175 Node v = _graph.target(e); 3176 int vb = _blossom_set>find(v); 3177 int vi = (*_node_index)[v]; 3178 3179 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot  3180 dualScale * _weight[e]; 3181 3182 if ((*_blossom_data)[vb].status == EVEN) { 3183 3184 int vt = _tree_set>find(vb); 3185 3186 typename std::map<int, Arc>::iterator it = 3187 (*_node_data)[ni].heap_index.find(vt); 3188 3189 if (it != (*_node_data)[ni].heap_index.end()) { 3190 if ((*_node_data)[ni].heap[it>second] > rw) { 3191 (*_node_data)[ni].heap.replace(it>second, e); 3192 (*_node_data)[ni].heap.decrease(e, rw); 3193 it>second = e; 3194 } 3195 } else { 3196 (*_node_data)[ni].heap.push(e, rw); 3197 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); 3198 } 3199 } 3200 } 3201 3202 if (!(*_node_data)[ni].heap.empty()) { 3203 _blossom_set>decrease(n, (*_node_data)[ni].heap.prio()); 3204 _delta2>push(nb, _blossom_set>classPrio(nb)); 3205 } 3206 } 3207 } 3208 3033 3209 /// \brief Start the algorithm 3034 3210 /// 3035 3211 /// This function starts the algorithm. 3036 3212 /// 3037 /// \pre \ref init() must be called before using this function. 3213 /// \pre \ref init() or \ref fractionalInit() must be called before 3214 /// using this function. 3038 3215 bool start() { 3039 3216 enum OpType { … … 3041 3218 }; 3042 3219 3043 int unmatched = _node_num; 3044 while (unmatched > 0) { 3220 if (_unmatched == 1) return false; 3221 3222 while (_unmatched > 0) { 3045 3223 Value d2 = !_delta2>empty() ? 3046 3224 _delta2>prio() : std::numeric_limits<Value>::max(); … … 3052 3230 _delta4>prio() : std::numeric_limits<Value>::max(); 3053 3231 3054 _delta_sum = d 2; OpType ot = D2;3055 if (d 3 < _delta_sum) { _delta_sum = d3; ot = D3; }3232 _delta_sum = d3; OpType ot = D3; 3233 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } 3056 3234 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } 3057 3235 … … 3086 3264 } else { 3087 3265 augmentOnEdge(e); 3088 unmatched = 2;3266 _unmatched = 2; 3089 3267 } 3090 3268 } … … 3105 3283 /// \note mwpm.run() is just a shortcut of the following code. 3106 3284 /// \code 3107 /// mwpm. init();3285 /// mwpm.fractionalInit(); 3108 3286 /// mwpm.start(); 3109 3287 /// \endcode 3110 3288 bool run() { 3111 init();3289 fractionalInit(); 3112 3290 return start(); 3113 3291 } … … 3116 3294 3117 3295 /// \name Primal Solution 3118 /// Functions to get the primal solution, i.e. the maximum weighted 3296 /// Functions to get the primal solution, i.e. the maximum weighted 3119 3297 /// perfect matching.\n 3120 3298 /// Either \ref run() or \ref start() function should be called before … … 3135 3313 } 3136 3314 } 3137 return sum / =2;3315 return sum / 2; 3138 3316 } 3139 3317 3140 3318 /// \brief Return \c true if the given edge is in the matching. 3141 3319 /// 3142 /// This function returns \c true if the given edge is in the found 3320 /// This function returns \c true if the given edge is in the found 3143 3321 /// matching. 3144 3322 /// … … 3151 3329 /// 3152 3330 /// This function returns the matching arc (or edge) incident to the 3153 /// given node in the found matching or \c INVALID if the node is 3331 /// given node in the found matching or \c INVALID if the node is 3154 3332 /// not covered by the matching. 3155 3333 /// … … 3169 3347 /// \brief Return the mate of the given node. 3170 3348 /// 3171 /// This function returns the mate of the given node in the found 3349 /// This function returns the mate of the given node in the found 3172 3350 /// matching or \c INVALID if the node is not covered by the matching. 3173 3351 /// … … 3188 3366 /// \brief Return the value of the dual solution. 3189 3367 /// 3190 /// This function returns the value of the dual solution. 3191 /// It should be equal to the primal value scaled by \ref dualScale 3368 /// This function returns the value of the dual solution. 3369 /// It should be equal to the primal value scaled by \ref dualScale 3192 3370 /// "dual scale". 3193 3371 /// … … 3244 3422 /// \brief Iterator for obtaining the nodes of a blossom. 3245 3423 /// 3246 /// This class provides an iterator for obtaining the nodes of the 3424 /// This class provides an iterator for obtaining the nodes of the 3247 3425 /// given blossom. It lists a subset of the nodes. 3248 /// Before using this iterator, you must allocate a 3426 /// Before using this iterator, you must allocate a 3249 3427 /// MaxWeightedPerfectMatching class and execute it. 3250 3428 class BlossomIt { … … 3255 3433 /// Constructor to get the nodes of the given variable. 3256 3434 /// 3257 /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 3258 /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 3435 /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 3436 /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 3259 3437 /// must be called before initializing this iterator. 3260 3438 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) … … 3302 3480 } //END OF NAMESPACE LEMON 3303 3481 3304 #endif //LEMON_MA X_MATCHING_H3482 #endif //LEMON_MATCHING_H 
lemon/matching.h
r949 r954 811 811 _matching = new MatchingMap(_graph); 812 812 } 813 813 814 if (!_node_potential) { 814 815 _node_potential = new NodePotential(_graph); 815 816 } 817 816 818 if (!_blossom_set) { 817 819 _blossom_index = new IntNodeMap(_graph); 818 820 _blossom_set = new BlossomSet(*_blossom_index); 819 821 _blossom_data = new RangeMap<BlossomData>(_blossom_num); 822 } else if (_blossom_data>size() != _blossom_num) { 823 delete _blossom_data; 824 _blossom_data = new RangeMap<BlossomData>(_blossom_num); 820 825 } 821 826 … … 824 829 _node_heap_index = new IntArcMap(_graph); 825 830 _node_data = new RangeMap<NodeData>(_node_num, 826 NodeData(*_node_heap_index)); 831 NodeData(*_node_heap_index)); 832 } else { 833 delete _node_data; 834 _node_data = new RangeMap<NodeData>(_node_num, 835 NodeData(*_node_heap_index)); 827 836 } 828 837 … … 830 839 _tree_set_index = new IntIntMap(_blossom_num); 831 840 _tree_set = new TreeSet(*_tree_set_index); 832 } 841 } else { 842 _tree_set_index>resize(_blossom_num); 843 } 844 833 845 if (!_delta1) { 834 846 _delta1_index = new IntNodeMap(_graph); 835 847 _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index); 836 848 } 849 837 850 if (!_delta2) { 838 851 _delta2_index = new IntIntMap(_blossom_num); 839 852 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index); 840 } 853 } else { 854 _delta2_index>resize(_blossom_num); 855 } 856 841 857 if (!_delta3) { 842 858 _delta3_index = new IntEdgeMap(_graph); 843 859 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index); 844 860 } 861 845 862 if (!_delta4) { 846 863 _delta4_index = new IntIntMap(_blossom_num); 847 864 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index); 865 } else { 866 _delta4_index>resize(_blossom_num); 848 867 } 849 868 } … … 1589 1608 createStructures(); 1590 1609 1610 _blossom_node_list.clear(); 1611 _blossom_potential.clear(); 1612 1591 1613 for (ArcIt e(_graph); e != INVALID; ++e) { 1592 1614 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; … … 1602 1624 (*_delta4_index)[i] = _delta4>PRE_HEAP; 1603 1625 } 1604 1626 1605 1627 _unmatched = _node_num; 1628 1629 _delta1>clear(); 1630 _delta2>clear(); 1631 _delta3>clear(); 1632 _delta4>clear(); 1633 _blossom_set>clear(); 1634 _tree_set>clear(); 1606 1635 1607 1636 int index = 0; … … 1615 1644 } 1616 1645 (*_node_index)[n] = index; 1646 (*_node_data)[index].heap_index.clear(); 1647 (*_node_data)[index].heap.clear(); 1617 1648 (*_node_data)[index].pot = max; 1618 1649 _delta1>push(n, max); … … 2238 2269 _matching = new MatchingMap(_graph); 2239 2270 } 2271 2240 2272 if (!_node_potential) { 2241 2273 _node_potential = new NodePotential(_graph); 2242 2274 } 2275 2243 2276 if (!_blossom_set) { 2244 2277 _blossom_index = new IntNodeMap(_graph); 2245 2278 _blossom_set = new BlossomSet(*_blossom_index); 2279 _blossom_data = new RangeMap<BlossomData>(_blossom_num); 2280 } else if (_blossom_data>size() != _blossom_num) { 2281 delete _blossom_data; 2246 2282 _blossom_data = new RangeMap<BlossomData>(_blossom_num); 2247 2283 } … … 2252 2288 _node_data = new RangeMap<NodeData>(_node_num, 2253 2289 NodeData(*_node_heap_index)); 2290 } else if (_node_data>size() != _node_num) { 2291 delete _node_data; 2292 _node_data = new RangeMap<NodeData>(_node_num, 2293 NodeData(*_node_heap_index)); 2254 2294 } 2255 2295 … … 2257 2297 _tree_set_index = new IntIntMap(_blossom_num); 2258 2298 _tree_set = new TreeSet(*_tree_set_index); 2259 } 2299 } else { 2300 _tree_set_index>resize(_blossom_num); 2301 } 2302 2260 2303 if (!_delta2) { 2261 2304 _delta2_index = new IntIntMap(_blossom_num); 2262 2305 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index); 2263 } 2306 } else { 2307 _delta2_index>resize(_blossom_num); 2308 } 2309 2264 2310 if (!_delta3) { 2265 2311 _delta3_index = new IntEdgeMap(_graph); 2266 2312 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index); 2267 2313 } 2314 2268 2315 if (!_delta4) { 2269 2316 _delta4_index = new IntIntMap(_blossom_num); 2270 2317 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index); 2318 } else { 2319 _delta4_index>resize(_blossom_num); 2271 2320 } 2272 2321 } … … 2969 3018 createStructures(); 2970 3019 3020 _blossom_node_list.clear(); 3021 _blossom_potential.clear(); 3022 2971 3023 for (ArcIt e(_graph); e != INVALID; ++e) { 2972 3024 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; … … 2981 3033 2982 3034 _unmatched = _node_num; 3035 3036 _delta2>clear(); 3037 _delta3>clear(); 3038 _delta4>clear(); 3039 _blossom_set>clear(); 3040 _tree_set>clear(); 2983 3041 2984 3042 int index = 0; … … 2992 3050 } 2993 3051 (*_node_index)[n] = index; 3052 (*_node_data)[index].heap_index.clear(); 3053 (*_node_data)[index].heap.clear(); 2994 3054 (*_node_data)[index].pot = max; 2995 3055 int blossom =
Note: See TracChangeset
for help on using the changeset viewer.