COIN-OR::LEMON - Graph Library

Changes in / [873:23da1b807280:874:d8ea85825e02] in lemon-1.2


Ignore:
Files:
2 added
6 edited

Legend:

Unmodified
Added
Removed
  • doc/groups.dox

    r873 r874  
    387387problem of maximum flow.
    388388
    389 \ref Circulation is a preflow push-relabel algorithm implemented directly 
     389\ref Circulation is a preflow push-relabel algorithm implemented directly
    390390for finding feasible circulations, which is a somewhat different problem,
    391391but it is strongly related to maximum flow.
     
    523523  Edmond's blossom shrinking algorithm for calculating maximum weighted
    524524  perfect matching in general graphs.
     525- \ref MaxFractionalMatching Push-relabel algorithm for calculating
     526  maximum cardinality fractional matching in general graphs.
     527- \ref MaxWeightedFractionalMatching Augmenting path algorithm for calculating
     528  maximum weighted fractional matching in general graphs.
     529- \ref MaxWeightedPerfectFractionalMatching
     530  Augmenting path algorithm for calculating maximum weighted
     531  perfect fractional matching in general graphs.
    525532
    526533\image html matching.png
  • lemon/Makefile.am

    r864 r874  
    8585        lemon/euler.h \
    8686        lemon/fib_heap.h \
     87        lemon/fractional_matching.h \
    8788        lemon/full_graph.h \
    8889        lemon/glpk.h \
  • lemon/matching.h

    r651 r870  
    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() {
     
    845850
    846851    void destroyStructures() {
    847       _node_num = countNodes(_graph);
    848       _blossom_num = _node_num * 3 / 2;
    849 
    850852      if (_matching) {
    851853        delete _matching;
     
    922924            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
    923925              _delta3->push(e, rw / 2);
    924             }
    925           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
    926             if (_delta3->state(e) != _delta3->IN_HEAP) {
    927               _delta3->push(e, rw);
    928926            }
    929927          } else {
     
    950948                               (*_blossom_data)[vb].offset);
    951949                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
    952                            (*_blossom_data)[vb].offset){
     950                           (*_blossom_data)[vb].offset) {
    953951                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
    954952                                   (*_blossom_data)[vb].offset);
     
    969967      if (!_blossom_set->trivial(blossom)) {
    970968        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
    971                      (*_blossom_data)[blossom].offset);
     969                      (*_blossom_data)[blossom].offset);
    972970      }
    973971    }
     
    10361034                }
    10371035              }
    1038             }
    1039 
    1040           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
    1041             if (_delta3->state(e) == _delta3->IN_HEAP) {
    1042               _delta3->erase(e);
    10431036            }
    10441037          } else {
     
    10781071          std::numeric_limits<Value>::max()) {
    10791072        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
    1080                        (*_blossom_data)[blossom].offset);
     1073                      (*_blossom_data)[blossom].offset);
    10811074      }
    10821075
     
    11171110            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
    11181111              _delta3->push(e, rw / 2);
    1119             }
    1120           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
    1121             if (_delta3->state(e) != _delta3->IN_HEAP) {
    1122               _delta3->push(e, rw);
    11231112            }
    11241113          } else {
     
    11581147    }
    11591148
    1160 
    1161     void matchedToUnmatched(int blossom) {
    1162       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
    1163         _delta2->erase(blossom);
    1164       }
    1165 
    1166       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
    1167            n != INVALID; ++n) {
    1168         int ni = (*_node_index)[n];
    1169 
    1170         _blossom_set->increase(n, std::numeric_limits<Value>::max());
    1171 
    1172         (*_node_data)[ni].heap.clear();
    1173         (*_node_data)[ni].heap_index.clear();
    1174 
    1175         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
    1176           Node v = _graph.target(e);
    1177           int vb = _blossom_set->find(v);
    1178           int vi = (*_node_index)[v];
    1179 
    1180           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
    1181             dualScale * _weight[e];
    1182 
    1183           if ((*_blossom_data)[vb].status == EVEN) {
    1184             if (_delta3->state(e) != _delta3->IN_HEAP) {
    1185               _delta3->push(e, rw);
    1186             }
    1187           }
    1188         }
    1189       }
    1190     }
    1191 
    1192     void unmatchedToMatched(int blossom) {
    1193       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
    1194            n != INVALID; ++n) {
    1195         int ni = (*_node_index)[n];
    1196 
    1197         for (InArcIt e(_graph, n); e != INVALID; ++e) {
    1198           Node v = _graph.source(e);
    1199           int vb = _blossom_set->find(v);
    1200           int vi = (*_node_index)[v];
    1201 
    1202           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
    1203             dualScale * _weight[e];
    1204 
    1205           if (vb == blossom) {
    1206             if (_delta3->state(e) == _delta3->IN_HEAP) {
    1207               _delta3->erase(e);
    1208             }
    1209           } else if ((*_blossom_data)[vb].status == EVEN) {
    1210 
    1211             if (_delta3->state(e) == _delta3->IN_HEAP) {
    1212               _delta3->erase(e);
    1213             }
    1214 
    1215             int vt = _tree_set->find(vb);
    1216 
    1217             Arc r = _graph.oppositeArc(e);
    1218 
    1219             typename std::map<int, Arc>::iterator it =
    1220               (*_node_data)[ni].heap_index.find(vt);
    1221 
    1222             if (it != (*_node_data)[ni].heap_index.end()) {
    1223               if ((*_node_data)[ni].heap[it->second] > rw) {
    1224                 (*_node_data)[ni].heap.replace(it->second, r);
    1225                 (*_node_data)[ni].heap.decrease(r, rw);
    1226                 it->second = r;
    1227               }
    1228             } else {
    1229               (*_node_data)[ni].heap.push(r, rw);
    1230               (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
    1231             }
    1232 
    1233             if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
    1234               _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
    1235 
    1236               if (_delta2->state(blossom) != _delta2->IN_HEAP) {
    1237                 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
    1238                              (*_blossom_data)[blossom].offset);
    1239               } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
    1240                          (*_blossom_data)[blossom].offset){
    1241                 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
    1242                                  (*_blossom_data)[blossom].offset);
    1243               }
    1244             }
    1245 
    1246           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
    1247             if (_delta3->state(e) == _delta3->IN_HEAP) {
    1248               _delta3->erase(e);
    1249             }
    1250           }
    1251         }
    1252       }
    1253     }
    1254 
    12551149    void alternatePath(int even, int tree) {
    12561150      int odd;
     
    12951189      destroyTree(tree);
    12961190
    1297       (*_blossom_data)[blossom].status = UNMATCHED;
    12981191      (*_blossom_data)[blossom].base = node;
    1299       matchedToUnmatched(blossom);
    1300     }
    1301 
     1192      (*_blossom_data)[blossom].next = INVALID;
     1193    }
    13021194
    13031195    void augmentOnEdge(const Edge& edge) {
     
    13061198      int right = _blossom_set->find(_graph.v(edge));
    13071199
    1308       if ((*_blossom_data)[left].status == EVEN) {
    1309         int left_tree = _tree_set->find(left);
    1310         alternatePath(left, left_tree);
    1311         destroyTree(left_tree);
    1312       } else {
    1313         (*_blossom_data)[left].status = MATCHED;
    1314         unmatchedToMatched(left);
    1315       }
    1316 
    1317       if ((*_blossom_data)[right].status == EVEN) {
    1318         int right_tree = _tree_set->find(right);
    1319         alternatePath(right, right_tree);
    1320         destroyTree(right_tree);
    1321       } else {
    1322         (*_blossom_data)[right].status = MATCHED;
    1323         unmatchedToMatched(right);
    1324       }
     1200      int left_tree = _tree_set->find(left);
     1201      alternatePath(left, left_tree);
     1202      destroyTree(left_tree);
     1203
     1204      int right_tree = _tree_set->find(right);
     1205      alternatePath(right, right_tree);
     1206      destroyTree(right_tree);
    13251207
    13261208      (*_blossom_data)[left].next = _graph.direct(edge, true);
    13271209      (*_blossom_data)[right].next = _graph.direct(edge, false);
     1210    }
     1211
     1212    void augmentOnArc(const Arc& arc) {
     1213
     1214      int left = _blossom_set->find(_graph.source(arc));
     1215      int right = _blossom_set->find(_graph.target(arc));
     1216
     1217      (*_blossom_data)[left].status = MATCHED;
     1218
     1219      int right_tree = _tree_set->find(right);
     1220      alternatePath(right, right_tree);
     1221      destroyTree(right_tree);
     1222
     1223      (*_blossom_data)[left].next = arc;
     1224      (*_blossom_data)[right].next = _graph.oppositeArc(arc);
    13281225    }
    13291226
     
    15301427          (*_blossom_data)[sb].pred = pred;
    15311428          (*_blossom_data)[sb].next =
    1532                            _graph.oppositeArc((*_blossom_data)[tb].next);
     1429            _graph.oppositeArc((*_blossom_data)[tb].next);
    15331430
    15341431          pred = (*_blossom_data)[ub].next;
     
    16301527
    16311528      for (int i = 0; i < int(blossoms.size()); ++i) {
    1632         if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
     1529        if ((*_blossom_data)[blossoms[i]].next != INVALID) {
    16331530
    16341531          Value offset = (*_blossom_data)[blossoms[i]].offset;
     
    16681565        _delta4_index(0), _delta4(0),
    16691566
    1670         _delta_sum() {}
     1567        _delta_sum(), _unmatched(0),
     1568
     1569        _fractional(0)
     1570    {}
    16711571
    16721572    ~MaxWeightedMatching() {
    16731573      destroyStructures();
     1574      if (_fractional) {
     1575        delete _fractional;
     1576      }
    16741577    }
    16751578
     
    16991602        (*_delta4_index)[i] = _delta4->PRE_HEAP;
    17001603      }
     1604
     1605      _unmatched = _node_num;
    17011606
    17021607      int index = 0;
     
    17341639    }
    17351640
     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
    17361778    /// \brief Start the algorithm
    17371779    ///
    17381780    /// This function starts the algorithm.
    17391781    ///
    1740     /// \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.
    17411784    void start() {
    17421785      enum OpType {
     
    17441787      };
    17451788
    1746       int unmatched = _node_num;
    1747       while (unmatched > 0) {
     1789      while (_unmatched > 0) {
    17481790        Value d1 = !_delta1->empty() ?
    17491791          _delta1->prio() : std::numeric_limits<Value>::max();
     
    17581800          _delta4->prio() : std::numeric_limits<Value>::max();
    17591801
    1760         _delta_sum = d1; OpType ot = D1;
     1802        _delta_sum = d3; OpType ot = D3;
     1803        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
    17611804        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
    1762         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
    17631805        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
    1764 
    17651806
    17661807        switch (ot) {
     
    17691810            Node n = _delta1->top();
    17701811            unmatchNode(n);
    1771             --unmatched;
     1812            --_unmatched;
    17721813          }
    17731814          break;
     
    17761817            int blossom = _delta2->top();
    17771818            Node n = _blossom_set->classTop(blossom);
    1778             Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
    1779             extendOnArc(e);
     1819            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
     1820            if ((*_blossom_data)[blossom].next == INVALID) {
     1821              augmentOnArc(a);
     1822              --_unmatched;
     1823            } else {
     1824              extendOnArc(a);
     1825            }
    17801826          }
    17811827          break;
     
    17901836              _delta3->pop();
    17911837            } else {
    1792               int left_tree;
    1793               if ((*_blossom_data)[left_blossom].status == EVEN) {
    1794                 left_tree = _tree_set->find(left_blossom);
    1795               } else {
    1796                 left_tree = -1;
    1797                 ++unmatched;
    1798               }
    1799               int right_tree;
    1800               if ((*_blossom_data)[right_blossom].status == EVEN) {
    1801                 right_tree = _tree_set->find(right_blossom);
    1802               } else {
    1803                 right_tree = -1;
    1804                 ++unmatched;
    1805               }
     1838              int left_tree = _tree_set->find(left_blossom);
     1839              int right_tree = _tree_set->find(right_blossom);
    18061840
    18071841              if (left_tree == right_tree) {
     
    18091843              } else {
    18101844                augmentOnEdge(e);
    1811                 unmatched -= 2;
     1845                _unmatched -= 2;
    18121846              }
    18131847            }
     
    18271861    /// \note mwm.run() is just a shortcut of the following code.
    18281862    /// \code
    1829     ///   mwm.init();
     1863    ///   mwm.fractionalInit();
    18301864    ///   mwm.start();
    18311865    /// \endcode
    18321866    void run() {
    1833       init();
     1867      fractionalInit();
    18341868      start();
    18351869    }
     
    18381872
    18391873    /// \name Primal Solution
    1840     /// Functions to get the primal solution, i.e. the maximum weighted 
     1874    /// Functions to get the primal solution, i.e. the maximum weighted
    18411875    /// matching.\n
    18421876    /// Either \ref run() or \ref start() function should be called before
     
    18571891        }
    18581892      }
    1859       return sum /= 2;
     1893      return sum / 2;
    18601894    }
    18611895
     
    18771911    /// \brief Return \c true if the given edge is in the matching.
    18781912    ///
    1879     /// This function returns \c true if the given edge is in the found 
     1913    /// This function returns \c true if the given edge is in the found
    18801914    /// matching.
    18811915    ///
     
    18881922    ///
    18891923    /// This function returns the matching arc (or edge) incident to the
    1890     /// given node in the found matching or \c INVALID if the node is 
     1924    /// given node in the found matching or \c INVALID if the node is
    18911925    /// not covered by the matching.
    18921926    ///
     
    19061940    /// \brief Return the mate of the given node.
    19071941    ///
    1908     /// This function returns the mate of the given node in the found 
     1942    /// This function returns the mate of the given node in the found
    19091943    /// matching or \c INVALID if the node is not covered by the matching.
    19101944    ///
     
    19261960    /// \brief Return the value of the dual solution.
    19271961    ///
    1928     /// This function returns the value of the dual solution. 
    1929     /// It should be equal to the primal value scaled by \ref dualScale 
     1962    /// This function returns the value of the dual solution.
     1963    /// It should be equal to the primal value scaled by \ref dualScale
    19301964    /// "dual scale".
    19311965    ///
     
    19822016    /// \brief Iterator for obtaining the nodes of a blossom.
    19832017    ///
    1984     /// This class provides an iterator for obtaining the nodes of the 
     2018    /// This class provides an iterator for obtaining the nodes of the
    19852019    /// given blossom. It lists a subset of the nodes.
    1986     /// Before using this iterator, you must allocate a 
     2020    /// Before using this iterator, you must allocate a
    19872021    /// MaxWeightedMatching class and execute it.
    19882022    class BlossomIt {
     
    19932027      /// Constructor to get the nodes of the given variable.
    19942028      ///
    1995       /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 
    1996       /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 
     2029      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
     2030      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
    19972031      /// called before initializing this iterator.
    19982032      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
     
    20472081  /// \f$O(nm\log n)\f$ time complexity.
    20482082  ///
    2049   /// The maximum weighted perfect matching problem is to find a subset of 
    2050   /// the edges in an undirected graph with maximum overall weight for which 
     2083  /// The maximum weighted perfect matching problem is to find a subset of
     2084  /// the edges in an undirected graph with maximum overall weight for which
    20512085  /// each node has exactly one incident edge.
    20522086  /// It can be formulated with the following linear program.
     
    20712105      \frac{\vert B \vert - 1}{2}z_B\f] */
    20722106  ///
    2073   /// The algorithm can be executed with the run() function. 
     2107  /// The algorithm can be executed with the run() function.
    20742108  /// After it the matching (the primal solution) and the dual solution
    2075   /// can be obtained using the query functions and the 
    2076   /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 
    2077   /// which is able to iterate on the nodes of a blossom. 
     2109  /// can be obtained using the query functions and the
     2110  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
     2111  /// which is able to iterate on the nodes of a blossom.
    20782112  /// If the value type is integer, then the dual solution is multiplied
    20792113  /// by \ref MaxWeightedMatching::dualScale "4".
    20802114  ///
    20812115  /// \tparam GR The undirected graph type the algorithm runs on.
    2082   /// \tparam WM The type edge weight map. The default type is 
     2116  /// \tparam WM The type edge weight map. The default type is
    20832117  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
    20842118#ifdef DOXYGEN
     
    21912225
    21922226    Value _delta_sum;
     2227    int _unmatched;
     2228
     2229    typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
     2230    FractionalMatching;
     2231    FractionalMatching *_fractional;
    21932232
    21942233    void createStructures() {
     
    22342273
    22352274    void destroyStructures() {
    2236       _node_num = countNodes(_graph);
    2237       _blossom_num = _node_num * 3 / 2;
    2238 
    22392275      if (_matching) {
    22402276        delete _matching;
     
    29092945        _delta4_index(0), _delta4(0),
    29102946
    2911         _delta_sum() {}
     2947        _delta_sum(), _unmatched(0),
     2948
     2949        _fractional(0)
     2950    {}
    29122951
    29132952    ~MaxWeightedPerfectMatching() {
    29142953      destroyStructures();
     2954      if (_fractional) {
     2955        delete _fractional;
     2956      }
    29152957    }
    29162958
     
    29372979        (*_delta4_index)[i] = _delta4->PRE_HEAP;
    29382980      }
     2981
     2982      _unmatched = _node_num;
    29392983
    29402984      int index = 0;
     
    29713015    }
    29723016
     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
    29733149    /// \brief Start the algorithm
    29743150    ///
    29753151    /// This function starts the algorithm.
    29763152    ///
    2977     /// \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.
    29783155    bool start() {
    29793156      enum OpType {
     
    29813158      };
    29823159
    2983       int unmatched = _node_num;
    2984       while (unmatched > 0) {
     3160      if (_unmatched == -1) return false;
     3161
     3162      while (_unmatched > 0) {
    29853163        Value d2 = !_delta2->empty() ?
    29863164          _delta2->prio() : std::numeric_limits<Value>::max();
     
    29923170          _delta4->prio() : std::numeric_limits<Value>::max();
    29933171
    2994         _delta_sum = d2; OpType ot = D2;
    2995         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
     3172        _delta_sum = d3; OpType ot = D3;
     3173        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
    29963174        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
    29973175
     
    30263204              } else {
    30273205                augmentOnEdge(e);
    3028                 unmatched -= 2;
     3206                _unmatched -= 2;
    30293207              }
    30303208            }
     
    30453223    /// \note mwpm.run() is just a shortcut of the following code.
    30463224    /// \code
    3047     ///   mwpm.init();
     3225    ///   mwpm.fractionalInit();
    30483226    ///   mwpm.start();
    30493227    /// \endcode
    30503228    bool run() {
    3051       init();
     3229      fractionalInit();
    30523230      return start();
    30533231    }
     
    30563234
    30573235    /// \name Primal Solution
    3058     /// Functions to get the primal solution, i.e. the maximum weighted 
     3236    /// Functions to get the primal solution, i.e. the maximum weighted
    30593237    /// perfect matching.\n
    30603238    /// Either \ref run() or \ref start() function should be called before
     
    30753253        }
    30763254      }
    3077       return sum /= 2;
     3255      return sum / 2;
    30783256    }
    30793257
    30803258    /// \brief Return \c true if the given edge is in the matching.
    30813259    ///
    3082     /// This function returns \c true if the given edge is in the found 
     3260    /// This function returns \c true if the given edge is in the found
    30833261    /// matching.
    30843262    ///
     
    30913269    ///
    30923270    /// This function returns the matching arc (or edge) incident to the
    3093     /// given node in the found matching or \c INVALID if the node is 
     3271    /// given node in the found matching or \c INVALID if the node is
    30943272    /// not covered by the matching.
    30953273    ///
     
    31093287    /// \brief Return the mate of the given node.
    31103288    ///
    3111     /// This function returns the mate of the given node in the found 
     3289    /// This function returns the mate of the given node in the found
    31123290    /// matching or \c INVALID if the node is not covered by the matching.
    31133291    ///
     
    31283306    /// \brief Return the value of the dual solution.
    31293307    ///
    3130     /// This function returns the value of the dual solution. 
    3131     /// It should be equal to the primal value scaled by \ref dualScale 
     3308    /// This function returns the value of the dual solution.
     3309    /// It should be equal to the primal value scaled by \ref dualScale
    31323310    /// "dual scale".
    31333311    ///
     
    31843362    /// \brief Iterator for obtaining the nodes of a blossom.
    31853363    ///
    3186     /// This class provides an iterator for obtaining the nodes of the 
     3364    /// This class provides an iterator for obtaining the nodes of the
    31873365    /// given blossom. It lists a subset of the nodes.
    3188     /// Before using this iterator, you must allocate a 
     3366    /// Before using this iterator, you must allocate a
    31893367    /// MaxWeightedPerfectMatching class and execute it.
    31903368    class BlossomIt {
     
    31953373      /// Constructor to get the nodes of the given variable.
    31963374      ///
    3197       /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 
    3198       /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 
     3375      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
     3376      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
    31993377      /// must be called before initializing this iterator.
    32003378      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
     
    32423420} //END OF NAMESPACE LEMON
    32433421
    3244 #endif //LEMON_MAX_MATCHING_H
     3422#endif //LEMON_MATCHING_H
  • test/CMakeLists.txt

    r799 r874  
    2222  error_test
    2323  euler_test
     24  fractional_matching_test
    2425  gomory_hu_test
    2526  graph_copy_test
  • test/Makefile.am

    r799 r874  
    2424        test/error_test \
    2525        test/euler_test \
     26        test/fractional_matching_test \
    2627        test/gomory_hu_test \
    2728        test/graph_copy_test \
     
    7273test_error_test_SOURCES = test/error_test.cc
    7374test_euler_test_SOURCES = test/euler_test.cc
     75test_fractional_matching_test_SOURCES = test/fractional_matching_test.cc
    7476test_gomory_hu_test_SOURCES = test/gomory_hu_test.cc
    7577test_graph_copy_test_SOURCES = test/graph_copy_test.cc
  • test/matching_test.cc

    r594 r870  
    402402      edgeMap("weight", weight).run();
    403403
    404     MaxMatching<SmartGraph> mm(graph);
    405     mm.run();
    406     checkMatching(graph, mm);
    407 
    408     MaxWeightedMatching<SmartGraph> mwm(graph, weight);
    409     mwm.run();
    410     checkWeightedMatching(graph, weight, mwm);
    411 
    412     MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
    413     bool perfect = mwpm.run();
    414 
    415     check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
    416           "Perfect matching found");
    417 
    418     if (perfect) {
    419       checkWeightedPerfectMatching(graph, weight, mwpm);
     404    bool perfect;
     405    {
     406      MaxMatching<SmartGraph> mm(graph);
     407      mm.run();
     408      checkMatching(graph, mm);
     409      perfect = 2 * mm.matchingSize() == countNodes(graph);
     410    }
     411
     412    {
     413      MaxWeightedMatching<SmartGraph> mwm(graph, weight);
     414      mwm.run();
     415      checkWeightedMatching(graph, weight, mwm);
     416    }
     417
     418    {
     419      MaxWeightedMatching<SmartGraph> mwm(graph, weight);
     420      mwm.init();
     421      mwm.start();
     422      checkWeightedMatching(graph, weight, mwm);
     423    }
     424
     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    }
     434
     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  }
Note: See TracChangeset for help on using the changeset viewer.