COIN-OR::LEMON - Graph Library

Changeset 875:07ec2b52e53d in lemon-1.2


Ignore:
Timestamp:
03/17/10 10:29:57 (9 years ago)
Author:
Alpar Juttner <alpar@…>
Branch:
default
Parents:
874:d8ea85825e02 (diff), 867:5b926cc36a4b (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Phase:
public
Message:

Merge #356

Location:
lemon
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • lemon/matching.h

    r867 r875  
    1717 */
    1818
    19 #ifndef LEMON_MAX_MATCHING_H
    20 #define LEMON_MAX_MATCHING_H
     19#ifndef LEMON_MATCHING_H
     20#define LEMON_MATCHING_H
    2121
    2222#include <vector>
     
    2929#include <lemon/bin_heap.h>
    3030#include <lemon/maps.h>
     31#include <lemon/fractional_matching.h>
    3132
    3233///\ingroup matching
     
    4243  /// This class implements Edmonds' alternating forest matching algorithm
    4344  /// 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
    4546  /// (the default is the empty one).
    4647  ///
     
    7071    ///\brief Status constants for Gallai-Edmonds decomposition.
    7172    ///
    72     ///These constants are used for indicating the Gallai-Edmonds 
     73    ///These constants are used for indicating the Gallai-Edmonds
    7374    ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
    7475    ///induce a subgraph with factor-critical components, the nodes with
    7576    ///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
    7778    ///perfect matching.
    7879    enum Status {
     
    513514    }
    514515
    515     /// \brief Start Edmonds' algorithm with a heuristic improvement 
     516    /// \brief Start Edmonds' algorithm with a heuristic improvement
    516517    /// for dense graphs
    517518    ///
     
    535536    /// \brief Run Edmonds' algorithm
    536537    ///
    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
    539540    /// (for which <tt>m>=2*n</tt> holds).
    540541    void run() {
     
    557558    /// \brief Return the size (cardinality) of the matching.
    558559    ///
    559     /// This function returns the size (cardinality) of the current matching. 
     560    /// This function returns the size (cardinality) of the current matching.
    560561    /// After run() it returns the size of the maximum matching in the graph.
    561562    int matchingSize() const {
     
    571572    /// \brief Return \c true if the given edge is in the matching.
    572573    ///
    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
    574575    /// matching.
    575576    bool matching(const Edge& edge) const {
     
    580581    ///
    581582    /// 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
    583584    /// not covered by the matching.
    584585    Arc matching(const Node& n) const {
     
    596597    /// \brief Return the mate of the given node.
    597598    ///
    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
    599600    /// matching or \c INVALID if the node is not covered by the matching.
    600601    Node mate(const Node& n) const {
     
    606607
    607608    /// \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
    609610    /// decomposition.
    610611
     
    649650  /// \f$O(nm\log n)\f$ time complexity.
    650651  ///
    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
    653654  /// each node has at most one incident edge.
    654655  /// It can be formulated with the following linear program.
     
    674675      \frac{\vert B \vert - 1}{2}z_B\f] */
    675676  ///
    676   /// The algorithm can be executed with the run() function. 
     677  /// The algorithm can be executed with the run() function.
    677678  /// 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.
    681682  /// If the value type is integer, then the dual solution is multiplied
    682683  /// by \ref MaxWeightedMatching::dualScale "4".
    683684  ///
    684685  /// \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
    686687  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
    687688#ifdef DOXYGEN
     
    746747
    747748    enum Status {
    748       EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
     749      EVEN = -1, MATCHED = 0, ODD = 1
    749750    };
    750751
     
    798799
    799800    Value _delta_sum;
     801    int _unmatched;
     802
     803    typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
     804    FractionalMatching *_fractional;
    800805
    801806    void createStructures() {
     
    864869
    865870    void destroyStructures() {
    866       _node_num = countNodes(_graph);
    867       _blossom_num = _node_num * 3 / 2;
    868 
    869871      if (_matching) {
    870872        delete _matching;
     
    941943            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
    942944              _delta3->push(e, rw / 2);
    943             }
    944           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
    945             if (_delta3->state(e) != _delta3->IN_HEAP) {
    946               _delta3->push(e, rw);
    947945            }
    948946          } else {
     
    969967                               (*_blossom_data)[vb].offset);
    970968                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
    971                            (*_blossom_data)[vb].offset){
     969                           (*_blossom_data)[vb].offset) {
    972970                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
    973971                                   (*_blossom_data)[vb].offset);
     
    988986      if (!_blossom_set->trivial(blossom)) {
    989987        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
    990                      (*_blossom_data)[blossom].offset);
     988                      (*_blossom_data)[blossom].offset);
    991989      }
    992990    }
     
    10551053                }
    10561054              }
    1057             }
    1058 
    1059           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
    1060             if (_delta3->state(e) == _delta3->IN_HEAP) {
    1061               _delta3->erase(e);
    10621055            }
    10631056          } else {
     
    10971090          std::numeric_limits<Value>::max()) {
    10981091        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
    1099                        (*_blossom_data)[blossom].offset);
     1092                      (*_blossom_data)[blossom].offset);
    11001093      }
    11011094
     
    11361129            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
    11371130              _delta3->push(e, rw / 2);
    1138             }
    1139           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
    1140             if (_delta3->state(e) != _delta3->IN_HEAP) {
    1141               _delta3->push(e, rw);
    11421131            }
    11431132          } else {
     
    11771166    }
    11781167
    1179 
    1180     void matchedToUnmatched(int blossom) {
    1181       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
    1182         _delta2->erase(blossom);
    1183       }
    1184 
    1185       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
    1186            n != INVALID; ++n) {
    1187         int ni = (*_node_index)[n];
    1188 
    1189         _blossom_set->increase(n, std::numeric_limits<Value>::max());
    1190 
    1191         (*_node_data)[ni].heap.clear();
    1192         (*_node_data)[ni].heap_index.clear();
    1193 
    1194         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
    1195           Node v = _graph.target(e);
    1196           int vb = _blossom_set->find(v);
    1197           int vi = (*_node_index)[v];
    1198 
    1199           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
    1200             dualScale * _weight[e];
    1201 
    1202           if ((*_blossom_data)[vb].status == EVEN) {
    1203             if (_delta3->state(e) != _delta3->IN_HEAP) {
    1204               _delta3->push(e, rw);
    1205             }
    1206           }
    1207         }
    1208       }
    1209     }
    1210 
    1211     void unmatchedToMatched(int blossom) {
    1212       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
    1213            n != INVALID; ++n) {
    1214         int ni = (*_node_index)[n];
    1215 
    1216         for (InArcIt e(_graph, n); e != INVALID; ++e) {
    1217           Node v = _graph.source(e);
    1218           int vb = _blossom_set->find(v);
    1219           int vi = (*_node_index)[v];
    1220 
    1221           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
    1222             dualScale * _weight[e];
    1223 
    1224           if (vb == blossom) {
    1225             if (_delta3->state(e) == _delta3->IN_HEAP) {
    1226               _delta3->erase(e);
    1227             }
    1228           } else if ((*_blossom_data)[vb].status == EVEN) {
    1229 
    1230             if (_delta3->state(e) == _delta3->IN_HEAP) {
    1231               _delta3->erase(e);
    1232             }
    1233 
    1234             int vt = _tree_set->find(vb);
    1235 
    1236             Arc r = _graph.oppositeArc(e);
    1237 
    1238             typename std::map<int, Arc>::iterator it =
    1239               (*_node_data)[ni].heap_index.find(vt);
    1240 
    1241             if (it != (*_node_data)[ni].heap_index.end()) {
    1242               if ((*_node_data)[ni].heap[it->second] > rw) {
    1243                 (*_node_data)[ni].heap.replace(it->second, r);
    1244                 (*_node_data)[ni].heap.decrease(r, rw);
    1245                 it->second = r;
    1246               }
    1247             } else {
    1248               (*_node_data)[ni].heap.push(r, rw);
    1249               (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
    1250             }
    1251 
    1252             if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
    1253               _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
    1254 
    1255               if (_delta2->state(blossom) != _delta2->IN_HEAP) {
    1256                 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
    1257                              (*_blossom_data)[blossom].offset);
    1258               } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
    1259                          (*_blossom_data)[blossom].offset){
    1260                 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
    1261                                  (*_blossom_data)[blossom].offset);
    1262               }
    1263             }
    1264 
    1265           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
    1266             if (_delta3->state(e) == _delta3->IN_HEAP) {
    1267               _delta3->erase(e);
    1268             }
    1269           }
    1270         }
    1271       }
    1272     }
    1273 
    12741168    void alternatePath(int even, int tree) {
    12751169      int odd;
     
    13141208      destroyTree(tree);
    13151209
    1316       (*_blossom_data)[blossom].status = UNMATCHED;
    13171210      (*_blossom_data)[blossom].base = node;
    1318       matchedToUnmatched(blossom);
    1319     }
    1320 
     1211      (*_blossom_data)[blossom].next = INVALID;
     1212    }
    13211213
    13221214    void augmentOnEdge(const Edge& edge) {
     
    13251217      int right = _blossom_set->find(_graph.v(edge));
    13261218
    1327       if ((*_blossom_data)[left].status == EVEN) {
    1328         int left_tree = _tree_set->find(left);
    1329         alternatePath(left, left_tree);
    1330         destroyTree(left_tree);
    1331       } else {
    1332         (*_blossom_data)[left].status = MATCHED;
    1333         unmatchedToMatched(left);
    1334       }
    1335 
    1336       if ((*_blossom_data)[right].status == EVEN) {
    1337         int right_tree = _tree_set->find(right);
    1338         alternatePath(right, right_tree);
    1339         destroyTree(right_tree);
    1340       } else {
    1341         (*_blossom_data)[right].status = MATCHED;
    1342         unmatchedToMatched(right);
    1343       }
     1219      int left_tree = _tree_set->find(left);
     1220      alternatePath(left, left_tree);
     1221      destroyTree(left_tree);
     1222
     1223      int right_tree = _tree_set->find(right);
     1224      alternatePath(right, right_tree);
     1225      destroyTree(right_tree);
    13441226
    13451227      (*_blossom_data)[left].next = _graph.direct(edge, true);
    13461228      (*_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);
    13471244    }
    13481245
     
    15491446          (*_blossom_data)[sb].pred = pred;
    15501447          (*_blossom_data)[sb].next =
    1551                            _graph.oppositeArc((*_blossom_data)[tb].next);
     1448            _graph.oppositeArc((*_blossom_data)[tb].next);
    15521449
    15531450          pred = (*_blossom_data)[ub].next;
     
    16491546
    16501547      for (int i = 0; i < int(blossoms.size()); ++i) {
    1651         if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
     1548        if ((*_blossom_data)[blossoms[i]].next != INVALID) {
    16521549
    16531550          Value offset = (*_blossom_data)[blossoms[i]].offset;
     
    16871584        _delta4_index(0), _delta4(0),
    16881585
    1689         _delta_sum() {}
     1586        _delta_sum(), _unmatched(0),
     1587
     1588        _fractional(0)
     1589    {}
    16901590
    16911591    ~MaxWeightedMatching() {
    16921592      destroyStructures();
     1593      if (_fractional) {
     1594        delete _fractional;
     1595      }
    16931596    }
    16941597
     
    17221625      }
    17231626     
     1627      _unmatched = _node_num;
     1628
    17241629      _delta1->clear();
    17251630      _delta2->clear();
     
    17651670    }
    17661671
     1672    /// \brief Initialize the algorithm with fractional matching
     1673    ///
     1674    /// This function initializes the algorithm with a fractional
     1675    /// matching. This initialization is also called jumpstart heuristic.
     1676    void fractionalInit() {
     1677      createStructures();
     1678     
     1679      if (_fractional == 0) {
     1680        _fractional = new FractionalMatching(_graph, _weight, false);
     1681      }
     1682      _fractional->run();
     1683
     1684      for (ArcIt e(_graph); e != INVALID; ++e) {
     1685        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
     1686      }
     1687      for (NodeIt n(_graph); n != INVALID; ++n) {
     1688        (*_delta1_index)[n] = _delta1->PRE_HEAP;
     1689      }
     1690      for (EdgeIt e(_graph); e != INVALID; ++e) {
     1691        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     1692      }
     1693      for (int i = 0; i < _blossom_num; ++i) {
     1694        (*_delta2_index)[i] = _delta2->PRE_HEAP;
     1695        (*_delta4_index)[i] = _delta4->PRE_HEAP;
     1696      }
     1697
     1698      _unmatched = 0;
     1699
     1700      int index = 0;
     1701      for (NodeIt n(_graph); n != INVALID; ++n) {
     1702        Value pot = _fractional->nodeValue(n);
     1703        (*_node_index)[n] = index;
     1704        (*_node_data)[index].pot = pot;
     1705        int blossom =
     1706          _blossom_set->insert(n, std::numeric_limits<Value>::max());
     1707
     1708        (*_blossom_data)[blossom].status = MATCHED;
     1709        (*_blossom_data)[blossom].pred = INVALID;
     1710        (*_blossom_data)[blossom].next = _fractional->matching(n);
     1711        if (_fractional->matching(n) == INVALID) {
     1712          (*_blossom_data)[blossom].base = n;
     1713        }
     1714        (*_blossom_data)[blossom].pot = 0;
     1715        (*_blossom_data)[blossom].offset = 0;
     1716        ++index;
     1717      }
     1718
     1719      typename Graph::template NodeMap<bool> processed(_graph, false);
     1720      for (NodeIt n(_graph); n != INVALID; ++n) {
     1721        if (processed[n]) continue;
     1722        processed[n] = true;
     1723        if (_fractional->matching(n) == INVALID) continue;
     1724        int num = 1;
     1725        Node v = _graph.target(_fractional->matching(n));
     1726        while (n != v) {
     1727          processed[v] = true;
     1728          v = _graph.target(_fractional->matching(v));
     1729          ++num;
     1730        }
     1731
     1732        if (num % 2 == 1) {
     1733          std::vector<int> subblossoms(num);
     1734
     1735          subblossoms[--num] = _blossom_set->find(n);
     1736          _delta1->push(n, _fractional->nodeValue(n));
     1737          v = _graph.target(_fractional->matching(n));
     1738          while (n != v) {
     1739            subblossoms[--num] = _blossom_set->find(v);
     1740            _delta1->push(v, _fractional->nodeValue(v));
     1741            v = _graph.target(_fractional->matching(v));           
     1742          }
     1743         
     1744          int surface =
     1745            _blossom_set->join(subblossoms.begin(), subblossoms.end());
     1746          (*_blossom_data)[surface].status = EVEN;
     1747          (*_blossom_data)[surface].pred = INVALID;
     1748          (*_blossom_data)[surface].next = INVALID;
     1749          (*_blossom_data)[surface].pot = 0;
     1750          (*_blossom_data)[surface].offset = 0;
     1751         
     1752          _tree_set->insert(surface);
     1753          ++_unmatched;
     1754        }
     1755      }
     1756
     1757      for (EdgeIt e(_graph); e != INVALID; ++e) {
     1758        int si = (*_node_index)[_graph.u(e)];
     1759        int sb = _blossom_set->find(_graph.u(e));
     1760        int ti = (*_node_index)[_graph.v(e)];
     1761        int tb = _blossom_set->find(_graph.v(e));
     1762        if ((*_blossom_data)[sb].status == EVEN &&
     1763            (*_blossom_data)[tb].status == EVEN && sb != tb) {
     1764          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
     1765                            dualScale * _weight[e]) / 2);
     1766        }
     1767      }
     1768
     1769      for (NodeIt n(_graph); n != INVALID; ++n) {
     1770        int nb = _blossom_set->find(n);
     1771        if ((*_blossom_data)[nb].status != MATCHED) continue;
     1772        int ni = (*_node_index)[n];
     1773
     1774        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
     1775          Node v = _graph.target(e);
     1776          int vb = _blossom_set->find(v);
     1777          int vi = (*_node_index)[v];
     1778
     1779          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
     1780            dualScale * _weight[e];
     1781
     1782          if ((*_blossom_data)[vb].status == EVEN) {
     1783
     1784            int vt = _tree_set->find(vb);
     1785
     1786            typename std::map<int, Arc>::iterator it =
     1787              (*_node_data)[ni].heap_index.find(vt);
     1788
     1789            if (it != (*_node_data)[ni].heap_index.end()) {
     1790              if ((*_node_data)[ni].heap[it->second] > rw) {
     1791                (*_node_data)[ni].heap.replace(it->second, e);
     1792                (*_node_data)[ni].heap.decrease(e, rw);
     1793                it->second = e;
     1794              }
     1795            } else {
     1796              (*_node_data)[ni].heap.push(e, rw);
     1797              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
     1798            }
     1799          }
     1800        }
     1801           
     1802        if (!(*_node_data)[ni].heap.empty()) {
     1803          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
     1804          _delta2->push(nb, _blossom_set->classPrio(nb));
     1805        }
     1806      }
     1807    }
     1808
    17671809    /// \brief Start the algorithm
    17681810    ///
    17691811    /// This function starts the algorithm.
    17701812    ///
    1771     /// \pre \ref init() must be called before using this function.
     1813    /// \pre \ref init() or \ref fractionalInit() must be called
     1814    /// before using this function.
    17721815    void start() {
    17731816      enum OpType {
     
    17751818      };
    17761819
    1777       int unmatched = _node_num;
    1778       while (unmatched > 0) {
     1820      while (_unmatched > 0) {
    17791821        Value d1 = !_delta1->empty() ?
    17801822          _delta1->prio() : std::numeric_limits<Value>::max();
     
    17891831          _delta4->prio() : std::numeric_limits<Value>::max();
    17901832
    1791         _delta_sum = d1; OpType ot = D1;
     1833        _delta_sum = d3; OpType ot = D3;
     1834        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
    17921835        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
    1793         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
    17941836        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
    1795 
    17961837
    17971838        switch (ot) {
     
    18001841            Node n = _delta1->top();
    18011842            unmatchNode(n);
    1802             --unmatched;
     1843            --_unmatched;
    18031844          }
    18041845          break;
     
    18071848            int blossom = _delta2->top();
    18081849            Node n = _blossom_set->classTop(blossom);
    1809             Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
    1810             extendOnArc(e);
     1850            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
     1851            if ((*_blossom_data)[blossom].next == INVALID) {
     1852              augmentOnArc(a);
     1853              --_unmatched;
     1854            } else {
     1855              extendOnArc(a);
     1856            }
    18111857          }
    18121858          break;
     
    18211867              _delta3->pop();
    18221868            } else {
    1823               int left_tree;
    1824               if ((*_blossom_data)[left_blossom].status == EVEN) {
    1825                 left_tree = _tree_set->find(left_blossom);
    1826               } else {
    1827                 left_tree = -1;
    1828                 ++unmatched;
    1829               }
    1830               int right_tree;
    1831               if ((*_blossom_data)[right_blossom].status == EVEN) {
    1832                 right_tree = _tree_set->find(right_blossom);
    1833               } else {
    1834                 right_tree = -1;
    1835                 ++unmatched;
    1836               }
     1869              int left_tree = _tree_set->find(left_blossom);
     1870              int right_tree = _tree_set->find(right_blossom);
    18371871
    18381872              if (left_tree == right_tree) {
     
    18401874              } else {
    18411875                augmentOnEdge(e);
    1842                 unmatched -= 2;
     1876                _unmatched -= 2;
    18431877              }
    18441878            }
     
    18581892    /// \note mwm.run() is just a shortcut of the following code.
    18591893    /// \code
    1860     ///   mwm.init();
     1894    ///   mwm.fractionalInit();
    18611895    ///   mwm.start();
    18621896    /// \endcode
    18631897    void run() {
    1864       init();
     1898      fractionalInit();
    18651899      start();
    18661900    }
     
    18691903
    18701904    /// \name Primal Solution
    1871     /// Functions to get the primal solution, i.e. the maximum weighted 
     1905    /// Functions to get the primal solution, i.e. the maximum weighted
    18721906    /// matching.\n
    18731907    /// Either \ref run() or \ref start() function should be called before
     
    18881922        }
    18891923      }
    1890       return sum /= 2;
     1924      return sum / 2;
    18911925    }
    18921926
     
    19081942    /// \brief Return \c true if the given edge is in the matching.
    19091943    ///
    1910     /// This function returns \c true if the given edge is in the found 
     1944    /// This function returns \c true if the given edge is in the found
    19111945    /// matching.
    19121946    ///
     
    19191953    ///
    19201954    /// This function returns the matching arc (or edge) incident to the
    1921     /// given node in the found matching or \c INVALID if the node is 
     1955    /// given node in the found matching or \c INVALID if the node is
    19221956    /// not covered by the matching.
    19231957    ///
     
    19371971    /// \brief Return the mate of the given node.
    19381972    ///
    1939     /// This function returns the mate of the given node in the found 
     1973    /// This function returns the mate of the given node in the found
    19401974    /// matching or \c INVALID if the node is not covered by the matching.
    19411975    ///
     
    19571991    /// \brief Return the value of the dual solution.
    19581992    ///
    1959     /// This function returns the value of the dual solution. 
    1960     /// It should be equal to the primal value scaled by \ref dualScale 
     1993    /// This function returns the value of the dual solution.
     1994    /// It should be equal to the primal value scaled by \ref dualScale
    19611995    /// "dual scale".
    19621996    ///
     
    20132047    /// \brief Iterator for obtaining the nodes of a blossom.
    20142048    ///
    2015     /// This class provides an iterator for obtaining the nodes of the 
     2049    /// This class provides an iterator for obtaining the nodes of the
    20162050    /// given blossom. It lists a subset of the nodes.
    2017     /// Before using this iterator, you must allocate a 
     2051    /// Before using this iterator, you must allocate a
    20182052    /// MaxWeightedMatching class and execute it.
    20192053    class BlossomIt {
     
    20242058      /// Constructor to get the nodes of the given variable.
    20252059      ///
    2026       /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 
    2027       /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 
     2060      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
     2061      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
    20282062      /// called before initializing this iterator.
    20292063      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
     
    20782112  /// \f$O(nm\log n)\f$ time complexity.
    20792113  ///
    2080   /// The maximum weighted perfect matching problem is to find a subset of 
    2081   /// the edges in an undirected graph with maximum overall weight for which 
     2114  /// The maximum weighted perfect matching problem is to find a subset of
     2115  /// the edges in an undirected graph with maximum overall weight for which
    20822116  /// each node has exactly one incident edge.
    20832117  /// It can be formulated with the following linear program.
     
    21022136      \frac{\vert B \vert - 1}{2}z_B\f] */
    21032137  ///
    2104   /// The algorithm can be executed with the run() function. 
     2138  /// The algorithm can be executed with the run() function.
    21052139  /// After it the matching (the primal solution) and the dual solution
    2106   /// can be obtained using the query functions and the 
    2107   /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 
    2108   /// which is able to iterate on the nodes of a blossom. 
     2140  /// can be obtained using the query functions and the
     2141  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
     2142  /// which is able to iterate on the nodes of a blossom.
    21092143  /// If the value type is integer, then the dual solution is multiplied
    21102144  /// by \ref MaxWeightedMatching::dualScale "4".
    21112145  ///
    21122146  /// \tparam GR The undirected graph type the algorithm runs on.
    2113   /// \tparam WM The type edge weight map. The default type is 
     2147  /// \tparam WM The type edge weight map. The default type is
    21142148  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
    21152149#ifdef DOXYGEN
     
    22222256
    22232257    Value _delta_sum;
     2258    int _unmatched;
     2259
     2260    typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
     2261    FractionalMatching;
     2262    FractionalMatching *_fractional;
    22242263
    22252264    void createStructures() {
     
    22832322
    22842323    void destroyStructures() {
    2285       _node_num = countNodes(_graph);
    2286       _blossom_num = _node_num * 3 / 2;
    2287 
    22882324      if (_matching) {
    22892325        delete _matching;
     
    29582994        _delta4_index(0), _delta4(0),
    29592995
    2960         _delta_sum() {}
     2996        _delta_sum(), _unmatched(0),
     2997
     2998        _fractional(0)
     2999    {}
    29613000
    29623001    ~MaxWeightedPerfectMatching() {
    29633002      destroyStructures();
     3003      if (_fractional) {
     3004        delete _fractional;
     3005      }
    29643006    }
    29653007
     
    29893031        (*_delta4_index)[i] = _delta4->PRE_HEAP;
    29903032      }
     3033
     3034      _unmatched = _node_num;
    29913035
    29923036      _delta2->clear();
     
    30313075    }
    30323076
     3077    /// \brief Initialize the algorithm with fractional matching
     3078    ///
     3079    /// This function initializes the algorithm with a fractional
     3080    /// matching. This initialization is also called jumpstart heuristic.
     3081    void fractionalInit() {
     3082      createStructures();
     3083     
     3084      if (_fractional == 0) {
     3085        _fractional = new FractionalMatching(_graph, _weight, false);
     3086      }
     3087      if (!_fractional->run()) {
     3088        _unmatched = -1;
     3089        return;
     3090      }
     3091
     3092      for (ArcIt e(_graph); e != INVALID; ++e) {
     3093        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
     3094      }
     3095      for (EdgeIt e(_graph); e != INVALID; ++e) {
     3096        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     3097      }
     3098      for (int i = 0; i < _blossom_num; ++i) {
     3099        (*_delta2_index)[i] = _delta2->PRE_HEAP;
     3100        (*_delta4_index)[i] = _delta4->PRE_HEAP;
     3101      }
     3102
     3103      _unmatched = 0;
     3104
     3105      int index = 0;
     3106      for (NodeIt n(_graph); n != INVALID; ++n) {
     3107        Value pot = _fractional->nodeValue(n);
     3108        (*_node_index)[n] = index;
     3109        (*_node_data)[index].pot = pot;
     3110        int blossom =
     3111          _blossom_set->insert(n, std::numeric_limits<Value>::max());
     3112
     3113        (*_blossom_data)[blossom].status = MATCHED;
     3114        (*_blossom_data)[blossom].pred = INVALID;
     3115        (*_blossom_data)[blossom].next = _fractional->matching(n);
     3116        (*_blossom_data)[blossom].pot = 0;
     3117        (*_blossom_data)[blossom].offset = 0;
     3118        ++index;
     3119      }
     3120
     3121      typename Graph::template NodeMap<bool> processed(_graph, false);
     3122      for (NodeIt n(_graph); n != INVALID; ++n) {
     3123        if (processed[n]) continue;
     3124        processed[n] = true;
     3125        if (_fractional->matching(n) == INVALID) continue;
     3126        int num = 1;
     3127        Node v = _graph.target(_fractional->matching(n));
     3128        while (n != v) {
     3129          processed[v] = true;
     3130          v = _graph.target(_fractional->matching(v));
     3131          ++num;
     3132        }
     3133
     3134        if (num % 2 == 1) {
     3135          std::vector<int> subblossoms(num);
     3136
     3137          subblossoms[--num] = _blossom_set->find(n);
     3138          v = _graph.target(_fractional->matching(n));
     3139          while (n != v) {
     3140            subblossoms[--num] = _blossom_set->find(v);
     3141            v = _graph.target(_fractional->matching(v));           
     3142          }
     3143         
     3144          int surface =
     3145            _blossom_set->join(subblossoms.begin(), subblossoms.end());
     3146          (*_blossom_data)[surface].status = EVEN;
     3147          (*_blossom_data)[surface].pred = INVALID;
     3148          (*_blossom_data)[surface].next = INVALID;
     3149          (*_blossom_data)[surface].pot = 0;
     3150          (*_blossom_data)[surface].offset = 0;
     3151         
     3152          _tree_set->insert(surface);
     3153          ++_unmatched;
     3154        }
     3155      }
     3156
     3157      for (EdgeIt e(_graph); e != INVALID; ++e) {
     3158        int si = (*_node_index)[_graph.u(e)];
     3159        int sb = _blossom_set->find(_graph.u(e));
     3160        int ti = (*_node_index)[_graph.v(e)];
     3161        int tb = _blossom_set->find(_graph.v(e));
     3162        if ((*_blossom_data)[sb].status == EVEN &&
     3163            (*_blossom_data)[tb].status == EVEN && sb != tb) {
     3164          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
     3165                            dualScale * _weight[e]) / 2);
     3166        }
     3167      }
     3168
     3169      for (NodeIt n(_graph); n != INVALID; ++n) {
     3170        int nb = _blossom_set->find(n);
     3171        if ((*_blossom_data)[nb].status != MATCHED) continue;
     3172        int ni = (*_node_index)[n];
     3173
     3174        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
     3175          Node v = _graph.target(e);
     3176          int vb = _blossom_set->find(v);
     3177          int vi = (*_node_index)[v];
     3178
     3179          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
     3180            dualScale * _weight[e];
     3181
     3182          if ((*_blossom_data)[vb].status == EVEN) {
     3183
     3184            int vt = _tree_set->find(vb);
     3185
     3186            typename std::map<int, Arc>::iterator it =
     3187              (*_node_data)[ni].heap_index.find(vt);
     3188
     3189            if (it != (*_node_data)[ni].heap_index.end()) {
     3190              if ((*_node_data)[ni].heap[it->second] > rw) {
     3191                (*_node_data)[ni].heap.replace(it->second, e);
     3192                (*_node_data)[ni].heap.decrease(e, rw);
     3193                it->second = e;
     3194              }
     3195            } else {
     3196              (*_node_data)[ni].heap.push(e, rw);
     3197              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
     3198            }
     3199          }
     3200        }
     3201           
     3202        if (!(*_node_data)[ni].heap.empty()) {
     3203          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
     3204          _delta2->push(nb, _blossom_set->classPrio(nb));
     3205        }
     3206      }
     3207    }
     3208
    30333209    /// \brief Start the algorithm
    30343210    ///
    30353211    /// This function starts the algorithm.
    30363212    ///
    3037     /// \pre \ref init() must be called before using this function.
     3213    /// \pre \ref init() or \ref fractionalInit() must be called before
     3214    /// using this function.
    30383215    bool start() {
    30393216      enum OpType {
     
    30413218      };
    30423219
    3043       int unmatched = _node_num;
    3044       while (unmatched > 0) {
     3220      if (_unmatched == -1) return false;
     3221
     3222      while (_unmatched > 0) {
    30453223        Value d2 = !_delta2->empty() ?
    30463224          _delta2->prio() : std::numeric_limits<Value>::max();
     
    30523230          _delta4->prio() : std::numeric_limits<Value>::max();
    30533231
    3054         _delta_sum = d2; OpType ot = D2;
    3055         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
     3232        _delta_sum = d3; OpType ot = D3;
     3233        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
    30563234        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
    30573235
     
    30863264              } else {
    30873265                augmentOnEdge(e);
    3088                 unmatched -= 2;
     3266                _unmatched -= 2;
    30893267              }
    30903268            }
     
    31053283    /// \note mwpm.run() is just a shortcut of the following code.
    31063284    /// \code
    3107     ///   mwpm.init();
     3285    ///   mwpm.fractionalInit();
    31083286    ///   mwpm.start();
    31093287    /// \endcode
    31103288    bool run() {
    3111       init();
     3289      fractionalInit();
    31123290      return start();
    31133291    }
     
    31163294
    31173295    /// \name Primal Solution
    3118     /// Functions to get the primal solution, i.e. the maximum weighted 
     3296    /// Functions to get the primal solution, i.e. the maximum weighted
    31193297    /// perfect matching.\n
    31203298    /// Either \ref run() or \ref start() function should be called before
     
    31353313        }
    31363314      }
    3137       return sum /= 2;
     3315      return sum / 2;
    31383316    }
    31393317
    31403318    /// \brief Return \c true if the given edge is in the matching.
    31413319    ///
    3142     /// This function returns \c true if the given edge is in the found 
     3320    /// This function returns \c true if the given edge is in the found
    31433321    /// matching.
    31443322    ///
     
    31513329    ///
    31523330    /// This function returns the matching arc (or edge) incident to the
    3153     /// given node in the found matching or \c INVALID if the node is 
     3331    /// given node in the found matching or \c INVALID if the node is
    31543332    /// not covered by the matching.
    31553333    ///
     
    31693347    /// \brief Return the mate of the given node.
    31703348    ///
    3171     /// This function returns the mate of the given node in the found 
     3349    /// This function returns the mate of the given node in the found
    31723350    /// matching or \c INVALID if the node is not covered by the matching.
    31733351    ///
     
    31883366    /// \brief Return the value of the dual solution.
    31893367    ///
    3190     /// This function returns the value of the dual solution. 
    3191     /// It should be equal to the primal value scaled by \ref dualScale 
     3368    /// This function returns the value of the dual solution.
     3369    /// It should be equal to the primal value scaled by \ref dualScale
    31923370    /// "dual scale".
    31933371    ///
     
    32443422    /// \brief Iterator for obtaining the nodes of a blossom.
    32453423    ///
    3246     /// This class provides an iterator for obtaining the nodes of the 
     3424    /// This class provides an iterator for obtaining the nodes of the
    32473425    /// given blossom. It lists a subset of the nodes.
    3248     /// Before using this iterator, you must allocate a 
     3426    /// Before using this iterator, you must allocate a
    32493427    /// MaxWeightedPerfectMatching class and execute it.
    32503428    class BlossomIt {
     
    32553433      /// Constructor to get the nodes of the given variable.
    32563434      ///
    3257       /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 
    3258       /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 
     3435      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
     3436      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
    32593437      /// must be called before initializing this iterator.
    32603438      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
     
    33023480} //END OF NAMESPACE LEMON
    33033481
    3304 #endif //LEMON_MAX_MATCHING_H
     3482#endif //LEMON_MATCHING_H
  • lemon/matching.h

    r870 r875  
    811811        _matching = new MatchingMap(_graph);
    812812      }
     813
    813814      if (!_node_potential) {
    814815        _node_potential = new NodePotential(_graph);
    815816      }
     817
    816818      if (!_blossom_set) {
    817819        _blossom_index = new IntNodeMap(_graph);
    818820        _blossom_set = new BlossomSet(*_blossom_index);
    819821        _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);
    820825      }
    821826
     
    824829        _node_heap_index = new IntArcMap(_graph);
    825830        _node_data = new RangeMap<NodeData>(_node_num,
    826                                               NodeData(*_node_heap_index));
     831                                            NodeData(*_node_heap_index));
     832      } else {
     833        delete _node_data;
     834        _node_data = new RangeMap<NodeData>(_node_num,
     835                                            NodeData(*_node_heap_index));
    827836      }
    828837
     
    830839        _tree_set_index = new IntIntMap(_blossom_num);
    831840        _tree_set = new TreeSet(*_tree_set_index);
    832       }
     841      } else {
     842        _tree_set_index->resize(_blossom_num);
     843      }
     844
    833845      if (!_delta1) {
    834846        _delta1_index = new IntNodeMap(_graph);
    835847        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
    836848      }
     849
    837850      if (!_delta2) {
    838851        _delta2_index = new IntIntMap(_blossom_num);
    839852        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
    840       }
     853      } else {
     854        _delta2_index->resize(_blossom_num);
     855      }
     856
    841857      if (!_delta3) {
    842858        _delta3_index = new IntEdgeMap(_graph);
    843859        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
    844860      }
     861
    845862      if (!_delta4) {
    846863        _delta4_index = new IntIntMap(_blossom_num);
    847864        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
     865      } else {
     866        _delta4_index->resize(_blossom_num);
    848867      }
    849868    }
     
    15891608      createStructures();
    15901609
     1610      _blossom_node_list.clear();
     1611      _blossom_potential.clear();
     1612
    15911613      for (ArcIt e(_graph); e != INVALID; ++e) {
    15921614        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
     
    16021624        (*_delta4_index)[i] = _delta4->PRE_HEAP;
    16031625      }
    1604 
     1626     
    16051627      _unmatched = _node_num;
     1628
     1629      _delta1->clear();
     1630      _delta2->clear();
     1631      _delta3->clear();
     1632      _delta4->clear();
     1633      _blossom_set->clear();
     1634      _tree_set->clear();
    16061635
    16071636      int index = 0;
     
    16151644        }
    16161645        (*_node_index)[n] = index;
     1646        (*_node_data)[index].heap_index.clear();
     1647        (*_node_data)[index].heap.clear();
    16171648        (*_node_data)[index].pot = max;
    16181649        _delta1->push(n, max);
     
    22382269        _matching = new MatchingMap(_graph);
    22392270      }
     2271
    22402272      if (!_node_potential) {
    22412273        _node_potential = new NodePotential(_graph);
    22422274      }
     2275
    22432276      if (!_blossom_set) {
    22442277        _blossom_index = new IntNodeMap(_graph);
    22452278        _blossom_set = new BlossomSet(*_blossom_index);
     2279        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
     2280      } else if (_blossom_data->size() != _blossom_num) {
     2281        delete _blossom_data;
    22462282        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
    22472283      }
     
    22522288        _node_data = new RangeMap<NodeData>(_node_num,
    22532289                                            NodeData(*_node_heap_index));
     2290      } else if (_node_data->size() != _node_num) {
     2291        delete _node_data;
     2292        _node_data = new RangeMap<NodeData>(_node_num,
     2293                                            NodeData(*_node_heap_index));
    22542294      }
    22552295
     
    22572297        _tree_set_index = new IntIntMap(_blossom_num);
    22582298        _tree_set = new TreeSet(*_tree_set_index);
    2259       }
     2299      } else {
     2300        _tree_set_index->resize(_blossom_num);
     2301      }
     2302
    22602303      if (!_delta2) {
    22612304        _delta2_index = new IntIntMap(_blossom_num);
    22622305        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
    2263       }
     2306      } else {
     2307        _delta2_index->resize(_blossom_num);
     2308      }
     2309
    22642310      if (!_delta3) {
    22652311        _delta3_index = new IntEdgeMap(_graph);
    22662312        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
    22672313      }
     2314
    22682315      if (!_delta4) {
    22692316        _delta4_index = new IntIntMap(_blossom_num);
    22702317        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
     2318      } else {
     2319        _delta4_index->resize(_blossom_num);
    22712320      }
    22722321    }
     
    29693018      createStructures();
    29703019
     3020      _blossom_node_list.clear();
     3021      _blossom_potential.clear();
     3022
    29713023      for (ArcIt e(_graph); e != INVALID; ++e) {
    29723024        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
     
    29813033
    29823034      _unmatched = _node_num;
     3035
     3036      _delta2->clear();
     3037      _delta3->clear();
     3038      _delta4->clear();
     3039      _blossom_set->clear();
     3040      _tree_set->clear();
    29833041
    29843042      int index = 0;
     
    29923050        }
    29933051        (*_node_index)[n] = index;
     3052        (*_node_data)[index].heap_index.clear();
     3053        (*_node_data)[index].heap.clear();
    29943054        (*_node_data)[index].pot = max;
    29953055        int blossom =
  • lemon/unionfind.h

    r867 r875  
    4444  /// This is a very simple but efficient implementation, providing
    4545  /// only four methods: join (union), find, insert and size.
    46   /// For more features see the \ref UnionFindEnum class.
     46  /// For more features, see the \ref UnionFindEnum class.
    4747  ///
    4848  /// It is primarily used in Kruskal algorithm for finding minimal
  • lemon/unionfind.h

    r800 r875  
    12891289        first_free_class(-1), first_free_node(-1) {}
    12901290
     1291    /// \brief Clears the union-find data structure
     1292    ///
     1293    /// Erase each item from the data structure.
     1294    void clear() {
     1295      nodes.clear();
     1296      classes.clear();
     1297      first_free_node = first_free_class = first_class = -1;
     1298    }
     1299
    12911300    /// \brief Insert a new node into a new component.
    12921301    ///
Note: See TracChangeset for help on using the changeset viewer.