Changeset 870:61120524af27 in lemon-1.2 for lemon
- Timestamp:
- 09/26/09 10:17:31 (15 years ago)
- Branch:
- default
- Phase:
- public
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
lemon/matching.h
r868 r870 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 … … 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() { … … 1560 1565 _delta4_index(0), _delta4(0), 1561 1566 1562 _delta_sum() {} 1567 _delta_sum(), _unmatched(0), 1568 1569 _fractional(0) 1570 {} 1563 1571 1564 1572 ~MaxWeightedMatching() { 1565 1573 destroyStructures(); 1574 if (_fractional) { 1575 delete _fractional; 1576 } 1566 1577 } 1567 1578 … … 1591 1602 (*_delta4_index)[i] = _delta4->PRE_HEAP; 1592 1603 } 1604 1605 _unmatched = _node_num; 1593 1606 1594 1607 int index = 0; … … 1626 1639 } 1627 1640 1641 /// \brief Initialize the algorithm with fractional matching 1642 /// 1643 /// This function initializes the algorithm with a fractional 1644 /// matching. This initialization is also called jumpstart heuristic. 1645 void fractionalInit() { 1646 createStructures(); 1647 1648 if (_fractional == 0) { 1649 _fractional = new FractionalMatching(_graph, _weight, false); 1650 } 1651 _fractional->run(); 1652 1653 for (ArcIt e(_graph); e != INVALID; ++e) { 1654 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; 1655 } 1656 for (NodeIt n(_graph); n != INVALID; ++n) { 1657 (*_delta1_index)[n] = _delta1->PRE_HEAP; 1658 } 1659 for (EdgeIt e(_graph); e != INVALID; ++e) { 1660 (*_delta3_index)[e] = _delta3->PRE_HEAP; 1661 } 1662 for (int i = 0; i < _blossom_num; ++i) { 1663 (*_delta2_index)[i] = _delta2->PRE_HEAP; 1664 (*_delta4_index)[i] = _delta4->PRE_HEAP; 1665 } 1666 1667 _unmatched = 0; 1668 1669 int index = 0; 1670 for (NodeIt n(_graph); n != INVALID; ++n) { 1671 Value pot = _fractional->nodeValue(n); 1672 (*_node_index)[n] = index; 1673 (*_node_data)[index].pot = pot; 1674 int blossom = 1675 _blossom_set->insert(n, std::numeric_limits<Value>::max()); 1676 1677 (*_blossom_data)[blossom].status = MATCHED; 1678 (*_blossom_data)[blossom].pred = INVALID; 1679 (*_blossom_data)[blossom].next = _fractional->matching(n); 1680 if (_fractional->matching(n) == INVALID) { 1681 (*_blossom_data)[blossom].base = n; 1682 } 1683 (*_blossom_data)[blossom].pot = 0; 1684 (*_blossom_data)[blossom].offset = 0; 1685 ++index; 1686 } 1687 1688 typename Graph::template NodeMap<bool> processed(_graph, false); 1689 for (NodeIt n(_graph); n != INVALID; ++n) { 1690 if (processed[n]) continue; 1691 processed[n] = true; 1692 if (_fractional->matching(n) == INVALID) continue; 1693 int num = 1; 1694 Node v = _graph.target(_fractional->matching(n)); 1695 while (n != v) { 1696 processed[v] = true; 1697 v = _graph.target(_fractional->matching(v)); 1698 ++num; 1699 } 1700 1701 if (num % 2 == 1) { 1702 std::vector<int> subblossoms(num); 1703 1704 subblossoms[--num] = _blossom_set->find(n); 1705 _delta1->push(n, _fractional->nodeValue(n)); 1706 v = _graph.target(_fractional->matching(n)); 1707 while (n != v) { 1708 subblossoms[--num] = _blossom_set->find(v); 1709 _delta1->push(v, _fractional->nodeValue(v)); 1710 v = _graph.target(_fractional->matching(v)); 1711 } 1712 1713 int surface = 1714 _blossom_set->join(subblossoms.begin(), subblossoms.end()); 1715 (*_blossom_data)[surface].status = EVEN; 1716 (*_blossom_data)[surface].pred = INVALID; 1717 (*_blossom_data)[surface].next = INVALID; 1718 (*_blossom_data)[surface].pot = 0; 1719 (*_blossom_data)[surface].offset = 0; 1720 1721 _tree_set->insert(surface); 1722 ++_unmatched; 1723 } 1724 } 1725 1726 for (EdgeIt e(_graph); e != INVALID; ++e) { 1727 int si = (*_node_index)[_graph.u(e)]; 1728 int sb = _blossom_set->find(_graph.u(e)); 1729 int ti = (*_node_index)[_graph.v(e)]; 1730 int tb = _blossom_set->find(_graph.v(e)); 1731 if ((*_blossom_data)[sb].status == EVEN && 1732 (*_blossom_data)[tb].status == EVEN && sb != tb) { 1733 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - 1734 dualScale * _weight[e]) / 2); 1735 } 1736 } 1737 1738 for (NodeIt n(_graph); n != INVALID; ++n) { 1739 int nb = _blossom_set->find(n); 1740 if ((*_blossom_data)[nb].status != MATCHED) continue; 1741 int ni = (*_node_index)[n]; 1742 1743 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 1744 Node v = _graph.target(e); 1745 int vb = _blossom_set->find(v); 1746 int vi = (*_node_index)[v]; 1747 1748 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 1749 dualScale * _weight[e]; 1750 1751 if ((*_blossom_data)[vb].status == EVEN) { 1752 1753 int vt = _tree_set->find(vb); 1754 1755 typename std::map<int, Arc>::iterator it = 1756 (*_node_data)[ni].heap_index.find(vt); 1757 1758 if (it != (*_node_data)[ni].heap_index.end()) { 1759 if ((*_node_data)[ni].heap[it->second] > rw) { 1760 (*_node_data)[ni].heap.replace(it->second, e); 1761 (*_node_data)[ni].heap.decrease(e, rw); 1762 it->second = e; 1763 } 1764 } else { 1765 (*_node_data)[ni].heap.push(e, rw); 1766 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); 1767 } 1768 } 1769 } 1770 1771 if (!(*_node_data)[ni].heap.empty()) { 1772 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); 1773 _delta2->push(nb, _blossom_set->classPrio(nb)); 1774 } 1775 } 1776 } 1777 1628 1778 /// \brief Start the algorithm 1629 1779 /// 1630 1780 /// This function starts the algorithm. 1631 1781 /// 1632 /// \pre \ref init() must be called before using this function. 1782 /// \pre \ref init() or \ref fractionalInit() must be called 1783 /// before using this function. 1633 1784 void start() { 1634 1785 enum OpType { … … 1636 1787 }; 1637 1788 1638 int unmatched = _node_num; 1639 while (unmatched > 0) { 1789 while (_unmatched > 0) { 1640 1790 Value d1 = !_delta1->empty() ? 1641 1791 _delta1->prio() : std::numeric_limits<Value>::max(); … … 1660 1810 Node n = _delta1->top(); 1661 1811 unmatchNode(n); 1662 -- unmatched;1812 --_unmatched; 1663 1813 } 1664 1814 break; … … 1670 1820 if ((*_blossom_data)[blossom].next == INVALID) { 1671 1821 augmentOnArc(a); 1672 -- unmatched;1822 --_unmatched; 1673 1823 } else { 1674 1824 extendOnArc(a); … … 1693 1843 } else { 1694 1844 augmentOnEdge(e); 1695 unmatched -= 2;1845 _unmatched -= 2; 1696 1846 } 1697 1847 } … … 1711 1861 /// \note mwm.run() is just a shortcut of the following code. 1712 1862 /// \code 1713 /// mwm. init();1863 /// mwm.fractionalInit(); 1714 1864 /// mwm.start(); 1715 1865 /// \endcode 1716 1866 void run() { 1717 init();1867 fractionalInit(); 1718 1868 start(); 1719 1869 } … … 2075 2225 2076 2226 Value _delta_sum; 2227 int _unmatched; 2228 2229 typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap> 2230 FractionalMatching; 2231 FractionalMatching *_fractional; 2077 2232 2078 2233 void createStructures() { … … 2790 2945 _delta4_index(0), _delta4(0), 2791 2946 2792 _delta_sum() {} 2947 _delta_sum(), _unmatched(0), 2948 2949 _fractional(0) 2950 {} 2793 2951 2794 2952 ~MaxWeightedPerfectMatching() { 2795 2953 destroyStructures(); 2954 if (_fractional) { 2955 delete _fractional; 2956 } 2796 2957 } 2797 2958 … … 2818 2979 (*_delta4_index)[i] = _delta4->PRE_HEAP; 2819 2980 } 2981 2982 _unmatched = _node_num; 2820 2983 2821 2984 int index = 0; … … 2852 3015 } 2853 3016 3017 /// \brief Initialize the algorithm with fractional matching 3018 /// 3019 /// This function initializes the algorithm with a fractional 3020 /// matching. This initialization is also called jumpstart heuristic. 3021 void fractionalInit() { 3022 createStructures(); 3023 3024 if (_fractional == 0) { 3025 _fractional = new FractionalMatching(_graph, _weight, false); 3026 } 3027 if (!_fractional->run()) { 3028 _unmatched = -1; 3029 return; 3030 } 3031 3032 for (ArcIt e(_graph); e != INVALID; ++e) { 3033 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; 3034 } 3035 for (EdgeIt e(_graph); e != INVALID; ++e) { 3036 (*_delta3_index)[e] = _delta3->PRE_HEAP; 3037 } 3038 for (int i = 0; i < _blossom_num; ++i) { 3039 (*_delta2_index)[i] = _delta2->PRE_HEAP; 3040 (*_delta4_index)[i] = _delta4->PRE_HEAP; 3041 } 3042 3043 _unmatched = 0; 3044 3045 int index = 0; 3046 for (NodeIt n(_graph); n != INVALID; ++n) { 3047 Value pot = _fractional->nodeValue(n); 3048 (*_node_index)[n] = index; 3049 (*_node_data)[index].pot = pot; 3050 int blossom = 3051 _blossom_set->insert(n, std::numeric_limits<Value>::max()); 3052 3053 (*_blossom_data)[blossom].status = MATCHED; 3054 (*_blossom_data)[blossom].pred = INVALID; 3055 (*_blossom_data)[blossom].next = _fractional->matching(n); 3056 (*_blossom_data)[blossom].pot = 0; 3057 (*_blossom_data)[blossom].offset = 0; 3058 ++index; 3059 } 3060 3061 typename Graph::template NodeMap<bool> processed(_graph, false); 3062 for (NodeIt n(_graph); n != INVALID; ++n) { 3063 if (processed[n]) continue; 3064 processed[n] = true; 3065 if (_fractional->matching(n) == INVALID) continue; 3066 int num = 1; 3067 Node v = _graph.target(_fractional->matching(n)); 3068 while (n != v) { 3069 processed[v] = true; 3070 v = _graph.target(_fractional->matching(v)); 3071 ++num; 3072 } 3073 3074 if (num % 2 == 1) { 3075 std::vector<int> subblossoms(num); 3076 3077 subblossoms[--num] = _blossom_set->find(n); 3078 v = _graph.target(_fractional->matching(n)); 3079 while (n != v) { 3080 subblossoms[--num] = _blossom_set->find(v); 3081 v = _graph.target(_fractional->matching(v)); 3082 } 3083 3084 int surface = 3085 _blossom_set->join(subblossoms.begin(), subblossoms.end()); 3086 (*_blossom_data)[surface].status = EVEN; 3087 (*_blossom_data)[surface].pred = INVALID; 3088 (*_blossom_data)[surface].next = INVALID; 3089 (*_blossom_data)[surface].pot = 0; 3090 (*_blossom_data)[surface].offset = 0; 3091 3092 _tree_set->insert(surface); 3093 ++_unmatched; 3094 } 3095 } 3096 3097 for (EdgeIt e(_graph); e != INVALID; ++e) { 3098 int si = (*_node_index)[_graph.u(e)]; 3099 int sb = _blossom_set->find(_graph.u(e)); 3100 int ti = (*_node_index)[_graph.v(e)]; 3101 int tb = _blossom_set->find(_graph.v(e)); 3102 if ((*_blossom_data)[sb].status == EVEN && 3103 (*_blossom_data)[tb].status == EVEN && sb != tb) { 3104 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - 3105 dualScale * _weight[e]) / 2); 3106 } 3107 } 3108 3109 for (NodeIt n(_graph); n != INVALID; ++n) { 3110 int nb = _blossom_set->find(n); 3111 if ((*_blossom_data)[nb].status != MATCHED) continue; 3112 int ni = (*_node_index)[n]; 3113 3114 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 3115 Node v = _graph.target(e); 3116 int vb = _blossom_set->find(v); 3117 int vi = (*_node_index)[v]; 3118 3119 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 3120 dualScale * _weight[e]; 3121 3122 if ((*_blossom_data)[vb].status == EVEN) { 3123 3124 int vt = _tree_set->find(vb); 3125 3126 typename std::map<int, Arc>::iterator it = 3127 (*_node_data)[ni].heap_index.find(vt); 3128 3129 if (it != (*_node_data)[ni].heap_index.end()) { 3130 if ((*_node_data)[ni].heap[it->second] > rw) { 3131 (*_node_data)[ni].heap.replace(it->second, e); 3132 (*_node_data)[ni].heap.decrease(e, rw); 3133 it->second = e; 3134 } 3135 } else { 3136 (*_node_data)[ni].heap.push(e, rw); 3137 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); 3138 } 3139 } 3140 } 3141 3142 if (!(*_node_data)[ni].heap.empty()) { 3143 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); 3144 _delta2->push(nb, _blossom_set->classPrio(nb)); 3145 } 3146 } 3147 } 3148 2854 3149 /// \brief Start the algorithm 2855 3150 /// 2856 3151 /// This function starts the algorithm. 2857 3152 /// 2858 /// \pre \ref init() must be called before using this function. 3153 /// \pre \ref init() or \ref fractionalInit() must be called before 3154 /// using this function. 2859 3155 bool start() { 2860 3156 enum OpType { … … 2862 3158 }; 2863 3159 2864 int unmatched = _node_num; 2865 while (unmatched > 0) { 3160 if (_unmatched == -1) return false; 3161 3162 while (_unmatched > 0) { 2866 3163 Value d2 = !_delta2->empty() ? 2867 3164 _delta2->prio() : std::numeric_limits<Value>::max(); … … 2907 3204 } else { 2908 3205 augmentOnEdge(e); 2909 unmatched -= 2;3206 _unmatched -= 2; 2910 3207 } 2911 3208 } … … 2926 3223 /// \note mwpm.run() is just a shortcut of the following code. 2927 3224 /// \code 2928 /// mwpm. init();3225 /// mwpm.fractionalInit(); 2929 3226 /// mwpm.start(); 2930 3227 /// \endcode 2931 3228 bool run() { 2932 init();3229 fractionalInit(); 2933 3230 return start(); 2934 3231 }
Note: See TracChangeset
for help on using the changeset viewer.