Changes in lemon/matching.h [651:3adf5e2d1e62:877:141f9c0db4a3] in lemon-1.2
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
lemon/matching.h
r651 r877 3 3 * This file is a part of LEMON, a generic C++ optimization library. 4 4 * 5 * Copyright (C) 2003-20 095 * Copyright (C) 2003-2010 6 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 7 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). … … 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 Gallai-Edmonds decomposition. 71 72 /// 72 ///These constants are used for indicating the Gallai-Edmonds 73 ///These constants are used for indicating the Gallai-Edmonds 73 74 ///decomposition of a graph. The nodes with status \c EVEN (or \c D) 74 75 ///induce a subgraph with factor-critical 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 Gallai-Edmonds 609 /// Functions to get the dual solution, i.e. the Gallai-Edmonds 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() { … … 806 811 _matching = new MatchingMap(_graph); 807 812 } 813 808 814 if (!_node_potential) { 809 815 _node_potential = new NodePotential(_graph); 810 816 } 817 811 818 if (!_blossom_set) { 812 819 _blossom_index = new IntNodeMap(_graph); 813 820 _blossom_set = new BlossomSet(*_blossom_index); 814 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); 815 825 } 816 826 … … 819 829 _node_heap_index = new IntArcMap(_graph); 820 830 _node_data = new RangeMap<NodeData>(_node_num, 821 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)); 822 836 } 823 837 … … 825 839 _tree_set_index = new IntIntMap(_blossom_num); 826 840 _tree_set = new TreeSet(*_tree_set_index); 827 } 841 } else { 842 _tree_set_index->resize(_blossom_num); 843 } 844 828 845 if (!_delta1) { 829 846 _delta1_index = new IntNodeMap(_graph); 830 847 _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index); 831 848 } 849 832 850 if (!_delta2) { 833 851 _delta2_index = new IntIntMap(_blossom_num); 834 852 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index); 835 } 853 } else { 854 _delta2_index->resize(_blossom_num); 855 } 856 836 857 if (!_delta3) { 837 858 _delta3_index = new IntEdgeMap(_graph); 838 859 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index); 839 860 } 861 840 862 if (!_delta4) { 841 863 _delta4_index = new IntIntMap(_blossom_num); 842 864 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index); 865 } else { 866 _delta4_index->resize(_blossom_num); 843 867 } 844 868 } 845 869 846 870 void destroyStructures() { 847 _node_num = countNodes(_graph);848 _blossom_num = _node_num * 3 / 2;849 850 871 if (_matching) { 851 872 delete _matching; … … 922 943 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { 923 944 _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 945 } 929 946 } else { … … 950 967 (*_blossom_data)[vb].offset); 951 968 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - 952 (*_blossom_data)[vb].offset) {969 (*_blossom_data)[vb].offset) { 953 970 _delta2->decrease(vb, _blossom_set->classPrio(vb) - 954 971 (*_blossom_data)[vb].offset); … … 969 986 if (!_blossom_set->trivial(blossom)) { 970 987 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + 971 (*_blossom_data)[blossom].offset);988 (*_blossom_data)[blossom].offset); 972 989 } 973 990 } … … 1036 1053 } 1037 1054 } 1038 }1039 1040 } else if ((*_blossom_data)[vb].status == UNMATCHED) {1041 if (_delta3->state(e) == _delta3->IN_HEAP) {1042 _delta3->erase(e);1043 1055 } 1044 1056 } else { … … 1078 1090 std::numeric_limits<Value>::max()) { 1079 1091 _delta2->push(blossom, _blossom_set->classPrio(blossom) - 1080 1092 (*_blossom_data)[blossom].offset); 1081 1093 } 1082 1094 … … 1117 1129 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { 1118 1130 _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 1131 } 1124 1132 } else { … … 1158 1166 } 1159 1167 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 1168 void alternatePath(int even, int tree) { 1256 1169 int odd; … … 1295 1208 destroyTree(tree); 1296 1209 1297 (*_blossom_data)[blossom].status = UNMATCHED;1298 1210 (*_blossom_data)[blossom].base = node; 1299 matchedToUnmatched(blossom); 1300 } 1301 1211 (*_blossom_data)[blossom].next = INVALID; 1212 } 1302 1213 1303 1214 void augmentOnEdge(const Edge& edge) { … … 1306 1217 int right = _blossom_set->find(_graph.v(edge)); 1307 1218 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 } 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); 1325 1226 1326 1227 (*_blossom_data)[left].next = _graph.direct(edge, true); 1327 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); 1328 1244 } 1329 1245 … … 1530 1446 (*_blossom_data)[sb].pred = pred; 1531 1447 (*_blossom_data)[sb].next = 1532 1448 _graph.oppositeArc((*_blossom_data)[tb].next); 1533 1449 1534 1450 pred = (*_blossom_data)[ub].next; … … 1630 1546 1631 1547 for (int i = 0; i < int(blossoms.size()); ++i) { 1632 if ((*_blossom_data)[blossoms[i]]. status == MATCHED) {1548 if ((*_blossom_data)[blossoms[i]].next != INVALID) { 1633 1549 1634 1550 Value offset = (*_blossom_data)[blossoms[i]].offset; … … 1668 1584 _delta4_index(0), _delta4(0), 1669 1585 1670 _delta_sum() {} 1586 _delta_sum(), _unmatched(0), 1587 1588 _fractional(0) 1589 {} 1671 1590 1672 1591 ~MaxWeightedMatching() { 1673 1592 destroyStructures(); 1593 if (_fractional) { 1594 delete _fractional; 1595 } 1674 1596 } 1675 1597 … … 1686 1608 createStructures(); 1687 1609 1610 _blossom_node_list.clear(); 1611 _blossom_potential.clear(); 1612 1688 1613 for (ArcIt e(_graph); e != INVALID; ++e) { 1689 1614 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; … … 1699 1624 (*_delta4_index)[i] = _delta4->PRE_HEAP; 1700 1625 } 1626 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(); 1701 1635 1702 1636 int index = 0; … … 1710 1644 } 1711 1645 (*_node_index)[n] = index; 1646 (*_node_data)[index].heap_index.clear(); 1647 (*_node_data)[index].heap.clear(); 1712 1648 (*_node_data)[index].pot = max; 1713 1649 _delta1->push(n, max); … … 1734 1670 } 1735 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 _blossom_node_list.clear(); 1680 _blossom_potential.clear(); 1681 1682 if (_fractional == 0) { 1683 _fractional = new FractionalMatching(_graph, _weight, false); 1684 } 1685 _fractional->run(); 1686 1687 for (ArcIt e(_graph); e != INVALID; ++e) { 1688 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; 1689 } 1690 for (NodeIt n(_graph); n != INVALID; ++n) { 1691 (*_delta1_index)[n] = _delta1->PRE_HEAP; 1692 } 1693 for (EdgeIt e(_graph); e != INVALID; ++e) { 1694 (*_delta3_index)[e] = _delta3->PRE_HEAP; 1695 } 1696 for (int i = 0; i < _blossom_num; ++i) { 1697 (*_delta2_index)[i] = _delta2->PRE_HEAP; 1698 (*_delta4_index)[i] = _delta4->PRE_HEAP; 1699 } 1700 1701 _unmatched = 0; 1702 1703 _delta1->clear(); 1704 _delta2->clear(); 1705 _delta3->clear(); 1706 _delta4->clear(); 1707 _blossom_set->clear(); 1708 _tree_set->clear(); 1709 1710 int index = 0; 1711 for (NodeIt n(_graph); n != INVALID; ++n) { 1712 Value pot = _fractional->nodeValue(n); 1713 (*_node_index)[n] = index; 1714 (*_node_data)[index].pot = pot; 1715 (*_node_data)[index].heap_index.clear(); 1716 (*_node_data)[index].heap.clear(); 1717 int blossom = 1718 _blossom_set->insert(n, std::numeric_limits<Value>::max()); 1719 1720 (*_blossom_data)[blossom].status = MATCHED; 1721 (*_blossom_data)[blossom].pred = INVALID; 1722 (*_blossom_data)[blossom].next = _fractional->matching(n); 1723 if (_fractional->matching(n) == INVALID) { 1724 (*_blossom_data)[blossom].base = n; 1725 } 1726 (*_blossom_data)[blossom].pot = 0; 1727 (*_blossom_data)[blossom].offset = 0; 1728 ++index; 1729 } 1730 1731 typename Graph::template NodeMap<bool> processed(_graph, false); 1732 for (NodeIt n(_graph); n != INVALID; ++n) { 1733 if (processed[n]) continue; 1734 processed[n] = true; 1735 if (_fractional->matching(n) == INVALID) continue; 1736 int num = 1; 1737 Node v = _graph.target(_fractional->matching(n)); 1738 while (n != v) { 1739 processed[v] = true; 1740 v = _graph.target(_fractional->matching(v)); 1741 ++num; 1742 } 1743 1744 if (num % 2 == 1) { 1745 std::vector<int> subblossoms(num); 1746 1747 subblossoms[--num] = _blossom_set->find(n); 1748 _delta1->push(n, _fractional->nodeValue(n)); 1749 v = _graph.target(_fractional->matching(n)); 1750 while (n != v) { 1751 subblossoms[--num] = _blossom_set->find(v); 1752 _delta1->push(v, _fractional->nodeValue(v)); 1753 v = _graph.target(_fractional->matching(v)); 1754 } 1755 1756 int surface = 1757 _blossom_set->join(subblossoms.begin(), subblossoms.end()); 1758 (*_blossom_data)[surface].status = EVEN; 1759 (*_blossom_data)[surface].pred = INVALID; 1760 (*_blossom_data)[surface].next = INVALID; 1761 (*_blossom_data)[surface].pot = 0; 1762 (*_blossom_data)[surface].offset = 0; 1763 1764 _tree_set->insert(surface); 1765 ++_unmatched; 1766 } 1767 } 1768 1769 for (EdgeIt e(_graph); e != INVALID; ++e) { 1770 int si = (*_node_index)[_graph.u(e)]; 1771 int sb = _blossom_set->find(_graph.u(e)); 1772 int ti = (*_node_index)[_graph.v(e)]; 1773 int tb = _blossom_set->find(_graph.v(e)); 1774 if ((*_blossom_data)[sb].status == EVEN && 1775 (*_blossom_data)[tb].status == EVEN && sb != tb) { 1776 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - 1777 dualScale * _weight[e]) / 2); 1778 } 1779 } 1780 1781 for (NodeIt n(_graph); n != INVALID; ++n) { 1782 int nb = _blossom_set->find(n); 1783 if ((*_blossom_data)[nb].status != MATCHED) continue; 1784 int ni = (*_node_index)[n]; 1785 1786 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 1787 Node v = _graph.target(e); 1788 int vb = _blossom_set->find(v); 1789 int vi = (*_node_index)[v]; 1790 1791 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 1792 dualScale * _weight[e]; 1793 1794 if ((*_blossom_data)[vb].status == EVEN) { 1795 1796 int vt = _tree_set->find(vb); 1797 1798 typename std::map<int, Arc>::iterator it = 1799 (*_node_data)[ni].heap_index.find(vt); 1800 1801 if (it != (*_node_data)[ni].heap_index.end()) { 1802 if ((*_node_data)[ni].heap[it->second] > rw) { 1803 (*_node_data)[ni].heap.replace(it->second, e); 1804 (*_node_data)[ni].heap.decrease(e, rw); 1805 it->second = e; 1806 } 1807 } else { 1808 (*_node_data)[ni].heap.push(e, rw); 1809 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); 1810 } 1811 } 1812 } 1813 1814 if (!(*_node_data)[ni].heap.empty()) { 1815 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); 1816 _delta2->push(nb, _blossom_set->classPrio(nb)); 1817 } 1818 } 1819 } 1820 1736 1821 /// \brief Start the algorithm 1737 1822 /// 1738 1823 /// This function starts the algorithm. 1739 1824 /// 1740 /// \pre \ref init() must be called before using this function. 1825 /// \pre \ref init() or \ref fractionalInit() must be called 1826 /// before using this function. 1741 1827 void start() { 1742 1828 enum OpType { … … 1744 1830 }; 1745 1831 1746 int unmatched = _node_num; 1747 while (unmatched > 0) { 1832 while (_unmatched > 0) { 1748 1833 Value d1 = !_delta1->empty() ? 1749 1834 _delta1->prio() : std::numeric_limits<Value>::max(); … … 1758 1843 _delta4->prio() : std::numeric_limits<Value>::max(); 1759 1844 1760 _delta_sum = d1; OpType ot = D1; 1845 _delta_sum = d3; OpType ot = D3; 1846 if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; } 1761 1847 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } 1762 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }1763 1848 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } 1764 1765 1849 1766 1850 switch (ot) { … … 1769 1853 Node n = _delta1->top(); 1770 1854 unmatchNode(n); 1771 -- unmatched;1855 --_unmatched; 1772 1856 } 1773 1857 break; … … 1776 1860 int blossom = _delta2->top(); 1777 1861 Node n = _blossom_set->classTop(blossom); 1778 Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); 1779 extendOnArc(e); 1862 Arc a = (*_node_data)[(*_node_index)[n]].heap.top(); 1863 if ((*_blossom_data)[blossom].next == INVALID) { 1864 augmentOnArc(a); 1865 --_unmatched; 1866 } else { 1867 extendOnArc(a); 1868 } 1780 1869 } 1781 1870 break; … … 1790 1879 _delta3->pop(); 1791 1880 } 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 } 1881 int left_tree = _tree_set->find(left_blossom); 1882 int right_tree = _tree_set->find(right_blossom); 1806 1883 1807 1884 if (left_tree == right_tree) { … … 1809 1886 } else { 1810 1887 augmentOnEdge(e); 1811 unmatched -= 2;1888 _unmatched -= 2; 1812 1889 } 1813 1890 } … … 1827 1904 /// \note mwm.run() is just a shortcut of the following code. 1828 1905 /// \code 1829 /// mwm. init();1906 /// mwm.fractionalInit(); 1830 1907 /// mwm.start(); 1831 1908 /// \endcode 1832 1909 void run() { 1833 init();1910 fractionalInit(); 1834 1911 start(); 1835 1912 } … … 1838 1915 1839 1916 /// \name Primal Solution 1840 /// Functions to get the primal solution, i.e. the maximum weighted 1917 /// Functions to get the primal solution, i.e. the maximum weighted 1841 1918 /// matching.\n 1842 1919 /// Either \ref run() or \ref start() function should be called before … … 1857 1934 } 1858 1935 } 1859 return sum / =2;1936 return sum / 2; 1860 1937 } 1861 1938 … … 1877 1954 /// \brief Return \c true if the given edge is in the matching. 1878 1955 /// 1879 /// This function returns \c true if the given edge is in the found 1956 /// This function returns \c true if the given edge is in the found 1880 1957 /// matching. 1881 1958 /// … … 1888 1965 /// 1889 1966 /// 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 1967 /// given node in the found matching or \c INVALID if the node is 1891 1968 /// not covered by the matching. 1892 1969 /// … … 1906 1983 /// \brief Return the mate of the given node. 1907 1984 /// 1908 /// This function returns the mate of the given node in the found 1985 /// This function returns the mate of the given node in the found 1909 1986 /// matching or \c INVALID if the node is not covered by the matching. 1910 1987 /// … … 1926 2003 /// \brief Return the value of the dual solution. 1927 2004 /// 1928 /// This function returns the value of the dual solution. 1929 /// It should be equal to the primal value scaled by \ref dualScale 2005 /// This function returns the value of the dual solution. 2006 /// It should be equal to the primal value scaled by \ref dualScale 1930 2007 /// "dual scale". 1931 2008 /// … … 1982 2059 /// \brief Iterator for obtaining the nodes of a blossom. 1983 2060 /// 1984 /// This class provides an iterator for obtaining the nodes of the 2061 /// This class provides an iterator for obtaining the nodes of the 1985 2062 /// given blossom. It lists a subset of the nodes. 1986 /// Before using this iterator, you must allocate a 2063 /// Before using this iterator, you must allocate a 1987 2064 /// MaxWeightedMatching class and execute it. 1988 2065 class BlossomIt { … … 1993 2070 /// Constructor to get the nodes of the given variable. 1994 2071 /// 1995 /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 1996 /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 2072 /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 2073 /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 1997 2074 /// called before initializing this iterator. 1998 2075 BlossomIt(const MaxWeightedMatching& algorithm, int variable) … … 2047 2124 /// \f$O(nm\log n)\f$ time complexity. 2048 2125 /// 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 2126 /// The maximum weighted perfect matching problem is to find a subset of 2127 /// the edges in an undirected graph with maximum overall weight for which 2051 2128 /// each node has exactly one incident edge. 2052 2129 /// It can be formulated with the following linear program. … … 2071 2148 \frac{\vert B \vert - 1}{2}z_B\f] */ 2072 2149 /// 2073 /// The algorithm can be executed with the run() function. 2150 /// The algorithm can be executed with the run() function. 2074 2151 /// 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. 2152 /// can be obtained using the query functions and the 2153 /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 2154 /// which is able to iterate on the nodes of a blossom. 2078 2155 /// If the value type is integer, then the dual solution is multiplied 2079 2156 /// by \ref MaxWeightedMatching::dualScale "4". 2080 2157 /// 2081 2158 /// \tparam GR The undirected graph type the algorithm runs on. 2082 /// \tparam WM The type edge weight map. The default type is 2159 /// \tparam WM The type edge weight map. The default type is 2083 2160 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>". 2084 2161 #ifdef DOXYGEN … … 2191 2268 2192 2269 Value _delta_sum; 2270 int _unmatched; 2271 2272 typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap> 2273 FractionalMatching; 2274 FractionalMatching *_fractional; 2193 2275 2194 2276 void createStructures() { … … 2199 2281 _matching = new MatchingMap(_graph); 2200 2282 } 2283 2201 2284 if (!_node_potential) { 2202 2285 _node_potential = new NodePotential(_graph); 2203 2286 } 2287 2204 2288 if (!_blossom_set) { 2205 2289 _blossom_index = new IntNodeMap(_graph); 2206 2290 _blossom_set = new BlossomSet(*_blossom_index); 2291 _blossom_data = new RangeMap<BlossomData>(_blossom_num); 2292 } else if (_blossom_data->size() != _blossom_num) { 2293 delete _blossom_data; 2207 2294 _blossom_data = new RangeMap<BlossomData>(_blossom_num); 2208 2295 } … … 2213 2300 _node_data = new RangeMap<NodeData>(_node_num, 2214 2301 NodeData(*_node_heap_index)); 2302 } else if (_node_data->size() != _node_num) { 2303 delete _node_data; 2304 _node_data = new RangeMap<NodeData>(_node_num, 2305 NodeData(*_node_heap_index)); 2215 2306 } 2216 2307 … … 2218 2309 _tree_set_index = new IntIntMap(_blossom_num); 2219 2310 _tree_set = new TreeSet(*_tree_set_index); 2220 } 2311 } else { 2312 _tree_set_index->resize(_blossom_num); 2313 } 2314 2221 2315 if (!_delta2) { 2222 2316 _delta2_index = new IntIntMap(_blossom_num); 2223 2317 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index); 2224 } 2318 } else { 2319 _delta2_index->resize(_blossom_num); 2320 } 2321 2225 2322 if (!_delta3) { 2226 2323 _delta3_index = new IntEdgeMap(_graph); 2227 2324 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index); 2228 2325 } 2326 2229 2327 if (!_delta4) { 2230 2328 _delta4_index = new IntIntMap(_blossom_num); 2231 2329 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index); 2330 } else { 2331 _delta4_index->resize(_blossom_num); 2232 2332 } 2233 2333 } 2234 2334 2235 2335 void destroyStructures() { 2236 _node_num = countNodes(_graph);2237 _blossom_num = _node_num * 3 / 2;2238 2239 2336 if (_matching) { 2240 2337 delete _matching; … … 2909 3006 _delta4_index(0), _delta4(0), 2910 3007 2911 _delta_sum() {} 3008 _delta_sum(), _unmatched(0), 3009 3010 _fractional(0) 3011 {} 2912 3012 2913 3013 ~MaxWeightedPerfectMatching() { 2914 3014 destroyStructures(); 3015 if (_fractional) { 3016 delete _fractional; 3017 } 2915 3018 } 2916 3019 … … 2927 3030 createStructures(); 2928 3031 3032 _blossom_node_list.clear(); 3033 _blossom_potential.clear(); 3034 2929 3035 for (ArcIt e(_graph); e != INVALID; ++e) { 2930 3036 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; … … 2937 3043 (*_delta4_index)[i] = _delta4->PRE_HEAP; 2938 3044 } 3045 3046 _unmatched = _node_num; 3047 3048 _delta2->clear(); 3049 _delta3->clear(); 3050 _delta4->clear(); 3051 _blossom_set->clear(); 3052 _tree_set->clear(); 2939 3053 2940 3054 int index = 0; … … 2948 3062 } 2949 3063 (*_node_index)[n] = index; 3064 (*_node_data)[index].heap_index.clear(); 3065 (*_node_data)[index].heap.clear(); 2950 3066 (*_node_data)[index].pot = max; 2951 3067 int blossom = … … 2971 3087 } 2972 3088 3089 /// \brief Initialize the algorithm with fractional matching 3090 /// 3091 /// This function initializes the algorithm with a fractional 3092 /// matching. This initialization is also called jumpstart heuristic. 3093 void fractionalInit() { 3094 createStructures(); 3095 3096 _blossom_node_list.clear(); 3097 _blossom_potential.clear(); 3098 3099 if (_fractional == 0) { 3100 _fractional = new FractionalMatching(_graph, _weight, false); 3101 } 3102 if (!_fractional->run()) { 3103 _unmatched = -1; 3104 return; 3105 } 3106 3107 for (ArcIt e(_graph); e != INVALID; ++e) { 3108 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; 3109 } 3110 for (EdgeIt e(_graph); e != INVALID; ++e) { 3111 (*_delta3_index)[e] = _delta3->PRE_HEAP; 3112 } 3113 for (int i = 0; i < _blossom_num; ++i) { 3114 (*_delta2_index)[i] = _delta2->PRE_HEAP; 3115 (*_delta4_index)[i] = _delta4->PRE_HEAP; 3116 } 3117 3118 _unmatched = 0; 3119 3120 _delta2->clear(); 3121 _delta3->clear(); 3122 _delta4->clear(); 3123 _blossom_set->clear(); 3124 _tree_set->clear(); 3125 3126 int index = 0; 3127 for (NodeIt n(_graph); n != INVALID; ++n) { 3128 Value pot = _fractional->nodeValue(n); 3129 (*_node_index)[n] = index; 3130 (*_node_data)[index].pot = pot; 3131 (*_node_data)[index].heap_index.clear(); 3132 (*_node_data)[index].heap.clear(); 3133 int blossom = 3134 _blossom_set->insert(n, std::numeric_limits<Value>::max()); 3135 3136 (*_blossom_data)[blossom].status = MATCHED; 3137 (*_blossom_data)[blossom].pred = INVALID; 3138 (*_blossom_data)[blossom].next = _fractional->matching(n); 3139 (*_blossom_data)[blossom].pot = 0; 3140 (*_blossom_data)[blossom].offset = 0; 3141 ++index; 3142 } 3143 3144 typename Graph::template NodeMap<bool> processed(_graph, false); 3145 for (NodeIt n(_graph); n != INVALID; ++n) { 3146 if (processed[n]) continue; 3147 processed[n] = true; 3148 if (_fractional->matching(n) == INVALID) continue; 3149 int num = 1; 3150 Node v = _graph.target(_fractional->matching(n)); 3151 while (n != v) { 3152 processed[v] = true; 3153 v = _graph.target(_fractional->matching(v)); 3154 ++num; 3155 } 3156 3157 if (num % 2 == 1) { 3158 std::vector<int> subblossoms(num); 3159 3160 subblossoms[--num] = _blossom_set->find(n); 3161 v = _graph.target(_fractional->matching(n)); 3162 while (n != v) { 3163 subblossoms[--num] = _blossom_set->find(v); 3164 v = _graph.target(_fractional->matching(v)); 3165 } 3166 3167 int surface = 3168 _blossom_set->join(subblossoms.begin(), subblossoms.end()); 3169 (*_blossom_data)[surface].status = EVEN; 3170 (*_blossom_data)[surface].pred = INVALID; 3171 (*_blossom_data)[surface].next = INVALID; 3172 (*_blossom_data)[surface].pot = 0; 3173 (*_blossom_data)[surface].offset = 0; 3174 3175 _tree_set->insert(surface); 3176 ++_unmatched; 3177 } 3178 } 3179 3180 for (EdgeIt e(_graph); e != INVALID; ++e) { 3181 int si = (*_node_index)[_graph.u(e)]; 3182 int sb = _blossom_set->find(_graph.u(e)); 3183 int ti = (*_node_index)[_graph.v(e)]; 3184 int tb = _blossom_set->find(_graph.v(e)); 3185 if ((*_blossom_data)[sb].status == EVEN && 3186 (*_blossom_data)[tb].status == EVEN && sb != tb) { 3187 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - 3188 dualScale * _weight[e]) / 2); 3189 } 3190 } 3191 3192 for (NodeIt n(_graph); n != INVALID; ++n) { 3193 int nb = _blossom_set->find(n); 3194 if ((*_blossom_data)[nb].status != MATCHED) continue; 3195 int ni = (*_node_index)[n]; 3196 3197 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 3198 Node v = _graph.target(e); 3199 int vb = _blossom_set->find(v); 3200 int vi = (*_node_index)[v]; 3201 3202 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 3203 dualScale * _weight[e]; 3204 3205 if ((*_blossom_data)[vb].status == EVEN) { 3206 3207 int vt = _tree_set->find(vb); 3208 3209 typename std::map<int, Arc>::iterator it = 3210 (*_node_data)[ni].heap_index.find(vt); 3211 3212 if (it != (*_node_data)[ni].heap_index.end()) { 3213 if ((*_node_data)[ni].heap[it->second] > rw) { 3214 (*_node_data)[ni].heap.replace(it->second, e); 3215 (*_node_data)[ni].heap.decrease(e, rw); 3216 it->second = e; 3217 } 3218 } else { 3219 (*_node_data)[ni].heap.push(e, rw); 3220 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); 3221 } 3222 } 3223 } 3224 3225 if (!(*_node_data)[ni].heap.empty()) { 3226 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); 3227 _delta2->push(nb, _blossom_set->classPrio(nb)); 3228 } 3229 } 3230 } 3231 2973 3232 /// \brief Start the algorithm 2974 3233 /// 2975 3234 /// This function starts the algorithm. 2976 3235 /// 2977 /// \pre \ref init() must be called before using this function. 3236 /// \pre \ref init() or \ref fractionalInit() must be called before 3237 /// using this function. 2978 3238 bool start() { 2979 3239 enum OpType { … … 2981 3241 }; 2982 3242 2983 int unmatched = _node_num; 2984 while (unmatched > 0) { 3243 if (_unmatched == -1) return false; 3244 3245 while (_unmatched > 0) { 2985 3246 Value d2 = !_delta2->empty() ? 2986 3247 _delta2->prio() : std::numeric_limits<Value>::max(); … … 2992 3253 _delta4->prio() : std::numeric_limits<Value>::max(); 2993 3254 2994 _delta_sum = d 2; OpType ot = D2;2995 if (d 3 < _delta_sum) { _delta_sum = d3; ot = D3; }3255 _delta_sum = d3; OpType ot = D3; 3256 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } 2996 3257 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } 2997 3258 … … 3026 3287 } else { 3027 3288 augmentOnEdge(e); 3028 unmatched -= 2;3289 _unmatched -= 2; 3029 3290 } 3030 3291 } … … 3045 3306 /// \note mwpm.run() is just a shortcut of the following code. 3046 3307 /// \code 3047 /// mwpm. init();3308 /// mwpm.fractionalInit(); 3048 3309 /// mwpm.start(); 3049 3310 /// \endcode 3050 3311 bool run() { 3051 init();3312 fractionalInit(); 3052 3313 return start(); 3053 3314 } … … 3056 3317 3057 3318 /// \name Primal Solution 3058 /// Functions to get the primal solution, i.e. the maximum weighted 3319 /// Functions to get the primal solution, i.e. the maximum weighted 3059 3320 /// perfect matching.\n 3060 3321 /// Either \ref run() or \ref start() function should be called before … … 3075 3336 } 3076 3337 } 3077 return sum / =2;3338 return sum / 2; 3078 3339 } 3079 3340 3080 3341 /// \brief Return \c true if the given edge is in the matching. 3081 3342 /// 3082 /// This function returns \c true if the given edge is in the found 3343 /// This function returns \c true if the given edge is in the found 3083 3344 /// matching. 3084 3345 /// … … 3091 3352 /// 3092 3353 /// 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 3354 /// given node in the found matching or \c INVALID if the node is 3094 3355 /// not covered by the matching. 3095 3356 /// … … 3109 3370 /// \brief Return the mate of the given node. 3110 3371 /// 3111 /// This function returns the mate of the given node in the found 3372 /// This function returns the mate of the given node in the found 3112 3373 /// matching or \c INVALID if the node is not covered by the matching. 3113 3374 /// … … 3128 3389 /// \brief Return the value of the dual solution. 3129 3390 /// 3130 /// This function returns the value of the dual solution. 3131 /// It should be equal to the primal value scaled by \ref dualScale 3391 /// This function returns the value of the dual solution. 3392 /// It should be equal to the primal value scaled by \ref dualScale 3132 3393 /// "dual scale". 3133 3394 /// … … 3184 3445 /// \brief Iterator for obtaining the nodes of a blossom. 3185 3446 /// 3186 /// This class provides an iterator for obtaining the nodes of the 3447 /// This class provides an iterator for obtaining the nodes of the 3187 3448 /// given blossom. It lists a subset of the nodes. 3188 /// Before using this iterator, you must allocate a 3449 /// Before using this iterator, you must allocate a 3189 3450 /// MaxWeightedPerfectMatching class and execute it. 3190 3451 class BlossomIt { … … 3195 3456 /// Constructor to get the nodes of the given variable. 3196 3457 /// 3197 /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 3198 /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 3458 /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 3459 /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 3199 3460 /// must be called before initializing this iterator. 3200 3461 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) … … 3242 3503 } //END OF NAMESPACE LEMON 3243 3504 3244 #endif //LEMON_MA X_MATCHING_H3505 #endif //LEMON_MATCHING_H
Note: See TracChangeset
for help on using the changeset viewer.