COIN-OR::LEMON - Graph Library

Ticket #314: 61120524af27.patch

File 61120524af27.patch, 16.4 KB (added by Balazs Dezso, 15 years ago)

Jumpstart heuristic for fractional matching

  • lemon/matching.h

    # HG changeset patch
    # User Balazs Dezso <deba@inf.elte.hu>
    # Date 1253953051 -7200
    # Node ID 61120524af27864803c669d537d0b9f2484abeb6
    # Parent  636dadefe1e6736f97988b05d6b4e7f95e4c5b9c
    Fractional matching initialization of weighted matchings (#314)
    
    diff -r 636dadefe1e6 -r 61120524af27 lemon/matching.h
    a b  
    2828#include <lemon/unionfind.h>
    2929#include <lemon/bin_heap.h>
    3030#include <lemon/maps.h>
     31#include <lemon/fractional_matching.h>
    3132
    3233///\ingroup matching
    3334///\file
     
    797798    BinHeap<Value, IntIntMap> *_delta4;
    798799
    799800    Value _delta_sum;
     801    int _unmatched;
     802
     803    typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
     804    FractionalMatching *_fractional;
    800805
    801806    void createStructures() {
    802807      _node_num = countNodes(_graph);
     
    15591564        _delta3_index(0), _delta3(0),
    15601565        _delta4_index(0), _delta4(0),
    15611566
    1562         _delta_sum() {}
     1567        _delta_sum(), _unmatched(0),
     1568
     1569        _fractional(0)
     1570    {}
    15631571
    15641572    ~MaxWeightedMatching() {
    15651573      destroyStructures();
     1574      if (_fractional) {
     1575        delete _fractional;
     1576      }
    15661577    }
    15671578
    15681579    /// \name Execution Control
     
    15911602        (*_delta4_index)[i] = _delta4->PRE_HEAP;
    15921603      }
    15931604
     1605      _unmatched = _node_num;
     1606
    15941607      int index = 0;
    15951608      for (NodeIt n(_graph); n != INVALID; ++n) {
    15961609        Value max = 0;
     
    16251638      }
    16261639    }
    16271640
     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
    16281778    /// \brief Start the algorithm
    16291779    ///
    16301780    /// This function starts the algorithm.
    16311781    ///
    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.
    16331784    void start() {
    16341785      enum OpType {
    16351786        D1, D2, D3, D4
    16361787      };
    16371788
    1638       int unmatched = _node_num;
    1639       while (unmatched > 0) {
     1789      while (_unmatched > 0) {
    16401790        Value d1 = !_delta1->empty() ?
    16411791          _delta1->prio() : std::numeric_limits<Value>::max();
    16421792
     
    16591809          {
    16601810            Node n = _delta1->top();
    16611811            unmatchNode(n);
    1662             --unmatched;
     1812            --_unmatched;
    16631813          }
    16641814          break;
    16651815        case D2:
     
    16691819            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
    16701820            if ((*_blossom_data)[blossom].next == INVALID) {
    16711821              augmentOnArc(a);
    1672               --unmatched;
     1822              --_unmatched;
    16731823            } else {
    16741824              extendOnArc(a);
    16751825            }
     
    16921842                shrinkOnEdge(e, left_tree);
    16931843              } else {
    16941844                augmentOnEdge(e);
    1695                 unmatched -= 2;
     1845                _unmatched -= 2;
    16961846              }
    16971847            }
    16981848          } break;
     
    17101860    ///
    17111861    /// \note mwm.run() is just a shortcut of the following code.
    17121862    /// \code
    1713     ///   mwm.init();
     1863    ///   mwm.fractionalInit();
    17141864    ///   mwm.start();
    17151865    /// \endcode
    17161866    void run() {
    1717       init();
     1867      fractionalInit();
    17181868      start();
    17191869    }
    17201870
     
    20742224    BinHeap<Value, IntIntMap> *_delta4;
    20752225
    20762226    Value _delta_sum;
     2227    int _unmatched;
     2228
     2229    typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
     2230    FractionalMatching;
     2231    FractionalMatching *_fractional;
    20772232
    20782233    void createStructures() {
    20792234      _node_num = countNodes(_graph);
     
    27892944        _delta3_index(0), _delta3(0),
    27902945        _delta4_index(0), _delta4(0),
    27912946
    2792         _delta_sum() {}
     2947        _delta_sum(), _unmatched(0),
     2948
     2949        _fractional(0)
     2950    {}
    27932951
    27942952    ~MaxWeightedPerfectMatching() {
    27952953      destroyStructures();
     2954      if (_fractional) {
     2955        delete _fractional;
     2956      }
    27962957    }
    27972958
    27982959    /// \name Execution Control
     
    28182979        (*_delta4_index)[i] = _delta4->PRE_HEAP;
    28192980      }
    28202981
     2982      _unmatched = _node_num;
     2983
    28212984      int index = 0;
    28222985      for (NodeIt n(_graph); n != INVALID; ++n) {
    28232986        Value max = - std::numeric_limits<Value>::max();
     
    28513014      }
    28523015    }
    28533016
     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
    28543149    /// \brief Start the algorithm
    28553150    ///
    28563151    /// This function starts the algorithm.
    28573152    ///
    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.
    28593155    bool start() {
    28603156      enum OpType {
    28613157        D2, D3, D4
    28623158      };
    28633159
    2864       int unmatched = _node_num;
    2865       while (unmatched > 0) {
     3160      if (_unmatched == -1) return false;
     3161
     3162      while (_unmatched > 0) {
    28663163        Value d2 = !_delta2->empty() ?
    28673164          _delta2->prio() : std::numeric_limits<Value>::max();
    28683165
     
    29063203                shrinkOnEdge(e, left_tree);
    29073204              } else {
    29083205                augmentOnEdge(e);
    2909                 unmatched -= 2;
     3206                _unmatched -= 2;
    29103207              }
    29113208            }
    29123209          } break;
     
    29253222    ///
    29263223    /// \note mwpm.run() is just a shortcut of the following code.
    29273224    /// \code
    2928     ///   mwpm.init();
     3225    ///   mwpm.fractionalInit();
    29293226    ///   mwpm.start();
    29303227    /// \endcode
    29313228    bool run() {
    2932       init();
     3229      fractionalInit();
    29333230      return start();
    29343231    }
    29353232
  • test/matching_test.cc

    diff -r 636dadefe1e6 -r 61120524af27 test/matching_test.cc
    a b  
    401401    graphReader(graph, lgfs).
    402402      edgeMap("weight", weight).run();
    403403
    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    }
    407411
    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    }
    411417
    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    }
    414424
    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    }
    417434
    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      }
    420444    }
    421445  }
    422446