gravatar
deba@inf.elte.hu
deba@inf.elte.hu
Fractional matching initialization of weighted matchings (#314)
0 2 0
default
2 files changed with 349 insertions and 28 deletions:
↑ Collapse diff ↑
Ignore white space 24 line context
... ...
@@ -19,24 +19,25 @@
19 19
#ifndef LEMON_MATCHING_H
20 20
#define LEMON_MATCHING_H
21 21

	
22 22
#include <vector>
23 23
#include <queue>
24 24
#include <set>
25 25
#include <limits>
26 26

	
27 27
#include <lemon/core.h>
28 28
#include <lemon/unionfind.h>
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
33 34
///\file
34 35
///\brief Maximum matching algorithms in general graphs.
35 36

	
36 37
namespace lemon {
37 38

	
38 39
  /// \ingroup matching
39 40
  ///
40 41
  /// \brief Maximum cardinality matching in general graphs
41 42
  ///
42 43
  /// This class implements Edmonds' alternating forest matching algorithm
... ...
@@ -788,24 +789,28 @@
788 789
    BinHeap<Value, IntNodeMap> *_delta1;
789 790

	
790 791
    IntIntMap *_delta2_index;
791 792
    BinHeap<Value, IntIntMap> *_delta2;
792 793

	
793 794
    IntEdgeMap *_delta3_index;
794 795
    BinHeap<Value, IntEdgeMap> *_delta3;
795 796

	
796 797
    IntIntMap *_delta4_index;
797 798
    BinHeap<Value, IntIntMap> *_delta4;
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() {
802 807
      _node_num = countNodes(_graph);
803 808
      _blossom_num = _node_num * 3 / 2;
804 809

	
805 810
      if (!_matching) {
806 811
        _matching = new MatchingMap(_graph);
807 812
      }
808 813
      if (!_node_potential) {
809 814
        _node_potential = new NodePotential(_graph);
810 815
      }
811 816
      if (!_blossom_set) {
... ...
@@ -1550,28 +1555,34 @@
1550 1555
        _node_potential(0), _blossom_potential(), _blossom_node_list(),
1551 1556
        _node_num(0), _blossom_num(0),
1552 1557

	
1553 1558
        _blossom_index(0), _blossom_set(0), _blossom_data(0),
1554 1559
        _node_index(0), _node_heap_index(0), _node_data(0),
1555 1560
        _tree_set_index(0), _tree_set(0),
1556 1561

	
1557 1562
        _delta1_index(0), _delta1(0),
1558 1563
        _delta2_index(0), _delta2(0),
1559 1564
        _delta3_index(0), _delta3(0),
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

	
1568 1579
    /// \name Execution Control
1569 1580
    /// The simplest way to execute the algorithm is to use the
1570 1581
    /// \ref run() member function.
1571 1582

	
1572 1583
    ///@{
1573 1584

	
1574 1585
    /// \brief Initialize the algorithm
1575 1586
    ///
1576 1587
    /// This function initializes the algorithm.
1577 1588
    void init() {
... ...
@@ -1582,24 +1593,26 @@
1582 1593
      }
1583 1594
      for (NodeIt n(_graph); n != INVALID; ++n) {
1584 1595
        (*_delta1_index)[n] = _delta1->PRE_HEAP;
1585 1596
      }
1586 1597
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1587 1598
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1588 1599
      }
1589 1600
      for (int i = 0; i < _blossom_num; ++i) {
1590 1601
        (*_delta2_index)[i] = _delta2->PRE_HEAP;
1591 1602
        (*_delta4_index)[i] = _delta4->PRE_HEAP;
1592 1603
      }
1593 1604

	
1605
      _unmatched = _node_num;
1606

	
1594 1607
      int index = 0;
1595 1608
      for (NodeIt n(_graph); n != INVALID; ++n) {
1596 1609
        Value max = 0;
1597 1610
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1598 1611
          if (_graph.target(e) == n) continue;
1599 1612
          if ((dualScale * _weight[e]) / 2 > max) {
1600 1613
            max = (dualScale * _weight[e]) / 2;
1601 1614
          }
1602 1615
        }
1603 1616
        (*_node_index)[n] = index;
1604 1617
        (*_node_data)[index].pot = max;
1605 1618
        _delta1->push(n, max);
... ...
@@ -1616,114 +1629,251 @@
1616 1629
        ++index;
1617 1630
      }
1618 1631
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1619 1632
        int si = (*_node_index)[_graph.u(e)];
1620 1633
        int ti = (*_node_index)[_graph.v(e)];
1621 1634
        if (_graph.u(e) != _graph.v(e)) {
1622 1635
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1623 1636
                            dualScale * _weight[e]) / 2);
1624 1637
        }
1625 1638
      }
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 {
1635 1786
        D1, D2, D3, D4
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();
1642 1792

	
1643 1793
        Value d2 = !_delta2->empty() ?
1644 1794
          _delta2->prio() : std::numeric_limits<Value>::max();
1645 1795

	
1646 1796
        Value d3 = !_delta3->empty() ?
1647 1797
          _delta3->prio() : std::numeric_limits<Value>::max();
1648 1798

	
1649 1799
        Value d4 = !_delta4->empty() ?
1650 1800
          _delta4->prio() : std::numeric_limits<Value>::max();
1651 1801

	
1652 1802
        _delta_sum = d3; OpType ot = D3;
1653 1803
        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1654 1804
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1655 1805
        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1656 1806

	
1657 1807
        switch (ot) {
1658 1808
        case D1:
1659 1809
          {
1660 1810
            Node n = _delta1->top();
1661 1811
            unmatchNode(n);
1662
            --unmatched;
1812
            --_unmatched;
1663 1813
          }
1664 1814
          break;
1665 1815
        case D2:
1666 1816
          {
1667 1817
            int blossom = _delta2->top();
1668 1818
            Node n = _blossom_set->classTop(blossom);
1669 1819
            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
1670 1820
            if ((*_blossom_data)[blossom].next == INVALID) {
1671 1821
              augmentOnArc(a);
1672
              --unmatched;
1822
              --_unmatched;
1673 1823
            } else {
1674 1824
              extendOnArc(a);
1675 1825
            }
1676 1826
          }
1677 1827
          break;
1678 1828
        case D3:
1679 1829
          {
1680 1830
            Edge e = _delta3->top();
1681 1831

	
1682 1832
            int left_blossom = _blossom_set->find(_graph.u(e));
1683 1833
            int right_blossom = _blossom_set->find(_graph.v(e));
1684 1834

	
1685 1835
            if (left_blossom == right_blossom) {
1686 1836
              _delta3->pop();
1687 1837
            } else {
1688 1838
              int left_tree = _tree_set->find(left_blossom);
1689 1839
              int right_tree = _tree_set->find(right_blossom);
1690 1840

	
1691 1841
              if (left_tree == right_tree) {
1692 1842
                shrinkOnEdge(e, left_tree);
1693 1843
              } else {
1694 1844
                augmentOnEdge(e);
1695
                unmatched -= 2;
1845
                _unmatched -= 2;
1696 1846
              }
1697 1847
            }
1698 1848
          } break;
1699 1849
        case D4:
1700 1850
          splitBlossom(_delta4->top());
1701 1851
          break;
1702 1852
        }
1703 1853
      }
1704 1854
      extractMatching();
1705 1855
    }
1706 1856

	
1707 1857
    /// \brief Run the algorithm.
1708 1858
    ///
1709 1859
    /// This method runs the \c %MaxWeightedMatching algorithm.
1710 1860
    ///
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
    }
1720 1870

	
1721 1871
    /// @}
1722 1872

	
1723 1873
    /// \name Primal Solution
1724 1874
    /// Functions to get the primal solution, i.e. the maximum weighted
1725 1875
    /// matching.\n
1726 1876
    /// Either \ref run() or \ref start() function should be called before
1727 1877
    /// using them.
1728 1878

	
1729 1879
    /// @{
... ...
@@ -2065,24 +2215,29 @@
2065 2215
    TreeSet *_tree_set;
2066 2216

	
2067 2217
    IntIntMap *_delta2_index;
2068 2218
    BinHeap<Value, IntIntMap> *_delta2;
2069 2219

	
2070 2220
    IntEdgeMap *_delta3_index;
2071 2221
    BinHeap<Value, IntEdgeMap> *_delta3;
2072 2222

	
2073 2223
    IntIntMap *_delta4_index;
2074 2224
    BinHeap<Value, IntIntMap> *_delta4;
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() {
2079 2234
      _node_num = countNodes(_graph);
2080 2235
      _blossom_num = _node_num * 3 / 2;
2081 2236

	
2082 2237
      if (!_matching) {
2083 2238
        _matching = new MatchingMap(_graph);
2084 2239
      }
2085 2240
      if (!_node_potential) {
2086 2241
        _node_potential = new NodePotential(_graph);
2087 2242
      }
2088 2243
      if (!_blossom_set) {
... ...
@@ -2780,28 +2935,34 @@
2780 2935
      : _graph(graph), _weight(weight), _matching(0),
2781 2936
        _node_potential(0), _blossom_potential(), _blossom_node_list(),
2782 2937
        _node_num(0), _blossom_num(0),
2783 2938

	
2784 2939
        _blossom_index(0), _blossom_set(0), _blossom_data(0),
2785 2940
        _node_index(0), _node_heap_index(0), _node_data(0),
2786 2941
        _tree_set_index(0), _tree_set(0),
2787 2942

	
2788 2943
        _delta2_index(0), _delta2(0),
2789 2944
        _delta3_index(0), _delta3(0),
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

	
2798 2959
    /// \name Execution Control
2799 2960
    /// The simplest way to execute the algorithm is to use the
2800 2961
    /// \ref run() member function.
2801 2962

	
2802 2963
    ///@{
2803 2964

	
2804 2965
    /// \brief Initialize the algorithm
2805 2966
    ///
2806 2967
    /// This function initializes the algorithm.
2807 2968
    void init() {
... ...
@@ -2809,24 +2970,26 @@
2809 2970

	
2810 2971
      for (ArcIt e(_graph); e != INVALID; ++e) {
2811 2972
        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
2812 2973
      }
2813 2974
      for (EdgeIt e(_graph); e != INVALID; ++e) {
2814 2975
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
2815 2976
      }
2816 2977
      for (int i = 0; i < _blossom_num; ++i) {
2817 2978
        (*_delta2_index)[i] = _delta2->PRE_HEAP;
2818 2979
        (*_delta4_index)[i] = _delta4->PRE_HEAP;
2819 2980
      }
2820 2981

	
2982
      _unmatched = _node_num;
2983

	
2821 2984
      int index = 0;
2822 2985
      for (NodeIt n(_graph); n != INVALID; ++n) {
2823 2986
        Value max = - std::numeric_limits<Value>::max();
2824 2987
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2825 2988
          if (_graph.target(e) == n) continue;
2826 2989
          if ((dualScale * _weight[e]) / 2 > max) {
2827 2990
            max = (dualScale * _weight[e]) / 2;
2828 2991
          }
2829 2992
        }
2830 2993
        (*_node_index)[n] = index;
2831 2994
        (*_node_data)[index].pot = max;
2832 2995
        int blossom =
... ...
@@ -2842,36 +3005,170 @@
2842 3005
        ++index;
2843 3006
      }
2844 3007
      for (EdgeIt e(_graph); e != INVALID; ++e) {
2845 3008
        int si = (*_node_index)[_graph.u(e)];
2846 3009
        int ti = (*_node_index)[_graph.v(e)];
2847 3010
        if (_graph.u(e) != _graph.v(e)) {
2848 3011
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
2849 3012
                            dualScale * _weight[e]) / 2);
2850 3013
        }
2851 3014
      }
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 {
2861 3157
        D2, D3, D4
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();
2868 3165

	
2869 3166
        Value d3 = !_delta3->empty() ?
2870 3167
          _delta3->prio() : std::numeric_limits<Value>::max();
2871 3168

	
2872 3169
        Value d4 = !_delta4->empty() ?
2873 3170
          _delta4->prio() : std::numeric_limits<Value>::max();
2874 3171

	
2875 3172
        _delta_sum = d3; OpType ot = D3;
2876 3173
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
2877 3174
        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
... ...
@@ -2897,48 +3194,48 @@
2897 3194
            int right_blossom = _blossom_set->find(_graph.v(e));
2898 3195

	
2899 3196
            if (left_blossom == right_blossom) {
2900 3197
              _delta3->pop();
2901 3198
            } else {
2902 3199
              int left_tree = _tree_set->find(left_blossom);
2903 3200
              int right_tree = _tree_set->find(right_blossom);
2904 3201

	
2905 3202
              if (left_tree == right_tree) {
2906 3203
                shrinkOnEdge(e, left_tree);
2907 3204
              } else {
2908 3205
                augmentOnEdge(e);
2909
                unmatched -= 2;
3206
                _unmatched -= 2;
2910 3207
              }
2911 3208
            }
2912 3209
          } break;
2913 3210
        case D4:
2914 3211
          splitBlossom(_delta4->top());
2915 3212
          break;
2916 3213
        }
2917 3214
      }
2918 3215
      extractMatching();
2919 3216
      return true;
2920 3217
    }
2921 3218

	
2922 3219
    /// \brief Run the algorithm.
2923 3220
    ///
2924 3221
    /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
2925 3222
    ///
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
    }
2935 3232

	
2936 3233
    /// @}
2937 3234

	
2938 3235
    /// \name Primal Solution
2939 3236
    /// Functions to get the primal solution, i.e. the maximum weighted
2940 3237
    /// perfect matching.\n
2941 3238
    /// Either \ref run() or \ref start() function should be called before
2942 3239
    /// using them.
2943 3240

	
2944 3241
    /// @{
Ignore white space 6 line context
... ...
@@ -392,33 +392,57 @@
392 392

	
393 393

	
394 394
int main() {
395 395

	
396 396
  for (int i = 0; i < lgfn; ++i) {
397 397
    SmartGraph graph;
398 398
    SmartGraph::EdgeMap<int> weight(graph);
399 399

	
400 400
    istringstream lgfs(lgf[i]);
401 401
    graphReader(graph, lgfs).
402 402
      edgeMap("weight", weight).run();
403 403

	
404
    MaxMatching<SmartGraph> mm(graph);
405
    mm.run();
406
    checkMatching(graph, mm);
404
    bool perfect;
405
    {
406
      MaxMatching<SmartGraph> mm(graph);
407
      mm.run();
408
      checkMatching(graph, mm);
409
      perfect = 2 * mm.matchingSize() == countNodes(graph);
410
    }
407 411

	
408
    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
409
    mwm.run();
410
    checkWeightedMatching(graph, weight, mwm);
412
    {
413
      MaxWeightedMatching<SmartGraph> mwm(graph, weight);
414
      mwm.run();
415
      checkWeightedMatching(graph, weight, mwm);
416
    }
411 417

	
412
    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
413
    bool perfect = mwpm.run();
418
    {
419
      MaxWeightedMatching<SmartGraph> mwm(graph, weight);
420
      mwm.init();
421
      mwm.start();
422
      checkWeightedMatching(graph, weight, mwm);
423
    }
414 424

	
415
    check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
416
          "Perfect matching found");
425
    {
426
      MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
427
      bool result = mwpm.run();
428
      
429
      check(result == perfect, "Perfect matching found");
430
      if (perfect) {
431
        checkWeightedPerfectMatching(graph, weight, mwpm);
432
      }
433
    }
417 434

	
418
    if (perfect) {
419
      checkWeightedPerfectMatching(graph, weight, mwpm);
435
    {
436
      MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
437
      mwpm.init();
438
      bool result = mwpm.start();
439
      
440
      check(result == perfect, "Perfect matching found");
441
      if (perfect) {
442
        checkWeightedPerfectMatching(graph, weight, mwpm);
443
      }
420 444
    }
421 445
  }
422 446

	
423 447
  return 0;
424 448
}
0 comments (0 inline)