COIN-OR::LEMON - Graph Library

Changeset 868:0513ccfea967 in lemon-main


Ignore:
Timestamp:
09/20/09 21:38:24 (15 years ago)
Author:
Balazs Dezso <deba@…>
Branch:
default
Phase:
public
Message:

General improvements in weighted matching algorithms (#314)

  • Fix include guard
  • Uniform handling of MATCHED and UNMATCHED blossoms
  • Prefer operations which decrease the number of trees
  • Fix improper use of '/='

The solved problems did not cause wrong solution.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/matching.h

    r651 r868  
    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>
     
    4242  /// This class implements Edmonds' alternating forest matching algorithm
    4343  /// for finding a maximum cardinality matching in a general undirected graph.
    44   /// It can be started from an arbitrary initial matching 
     44  /// It can be started from an arbitrary initial matching
    4545  /// (the default is the empty one).
    4646  ///
     
    7070    ///\brief Status constants for Gallai-Edmonds decomposition.
    7171    ///
    72     ///These constants are used for indicating the Gallai-Edmonds 
     72    ///These constants are used for indicating the Gallai-Edmonds
    7373    ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
    7474    ///induce a subgraph with factor-critical components, the nodes with
    7575    ///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 
     76    ///with status \c MATCHED (or \c C) induce a subgraph having a
    7777    ///perfect matching.
    7878    enum Status {
     
    513513    }
    514514
    515     /// \brief Start Edmonds' algorithm with a heuristic improvement 
     515    /// \brief Start Edmonds' algorithm with a heuristic improvement
    516516    /// for dense graphs
    517517    ///
     
    535535    /// \brief Run Edmonds' algorithm
    536536    ///
    537     /// This function runs Edmonds' algorithm. An additional heuristic of 
    538     /// postponing shrinks is used for relatively dense graphs 
     537    /// This function runs Edmonds' algorithm. An additional heuristic of
     538    /// postponing shrinks is used for relatively dense graphs
    539539    /// (for which <tt>m>=2*n</tt> holds).
    540540    void run() {
     
    557557    /// \brief Return the size (cardinality) of the matching.
    558558    ///
    559     /// This function returns the size (cardinality) of the current matching. 
     559    /// This function returns the size (cardinality) of the current matching.
    560560    /// After run() it returns the size of the maximum matching in the graph.
    561561    int matchingSize() const {
     
    571571    /// \brief Return \c true if the given edge is in the matching.
    572572    ///
    573     /// This function returns \c true if the given edge is in the current 
     573    /// This function returns \c true if the given edge is in the current
    574574    /// matching.
    575575    bool matching(const Edge& edge) const {
     
    580580    ///
    581581    /// 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 
     582    /// given node in the current matching or \c INVALID if the node is
    583583    /// not covered by the matching.
    584584    Arc matching(const Node& n) const {
     
    596596    /// \brief Return the mate of the given node.
    597597    ///
    598     /// This function returns the mate of the given node in the current 
     598    /// This function returns the mate of the given node in the current
    599599    /// matching or \c INVALID if the node is not covered by the matching.
    600600    Node mate(const Node& n) const {
     
    606606
    607607    /// \name Dual Solution
    608     /// Functions to get the dual solution, i.e. the Gallai-Edmonds 
     608    /// Functions to get the dual solution, i.e. the Gallai-Edmonds
    609609    /// decomposition.
    610610
     
    649649  /// \f$O(nm\log n)\f$ time complexity.
    650650  ///
    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 
     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
    653653  /// each node has at most one incident edge.
    654654  /// It can be formulated with the following linear program.
     
    674674      \frac{\vert B \vert - 1}{2}z_B\f] */
    675675  ///
    676   /// The algorithm can be executed with the run() function. 
     676  /// The algorithm can be executed with the run() function.
    677677  /// 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. 
     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.
    681681  /// If the value type is integer, then the dual solution is multiplied
    682682  /// by \ref MaxWeightedMatching::dualScale "4".
    683683  ///
    684684  /// \tparam GR The undirected graph type the algorithm runs on.
    685   /// \tparam WM The type edge weight map. The default type is 
     685  /// \tparam WM The type edge weight map. The default type is
    686686  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
    687687#ifdef DOXYGEN
     
    746746
    747747    enum Status {
    748       EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
     748      EVEN = -1, MATCHED = 0, ODD = 1
    749749    };
    750750
     
    845845
    846846    void destroyStructures() {
    847       _node_num = countNodes(_graph);
    848       _blossom_num = _node_num * 3 / 2;
    849 
    850847      if (_matching) {
    851848        delete _matching;
     
    922919            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
    923920              _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);
    928921            }
    929922          } else {
     
    950943                               (*_blossom_data)[vb].offset);
    951944                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
    952                            (*_blossom_data)[vb].offset){
     945                           (*_blossom_data)[vb].offset) {
    953946                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
    954947                                   (*_blossom_data)[vb].offset);
     
    969962      if (!_blossom_set->trivial(blossom)) {
    970963        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
    971                      (*_blossom_data)[blossom].offset);
     964                      (*_blossom_data)[blossom].offset);
    972965      }
    973966    }
     
    10371030              }
    10381031            }
    1039 
    1040           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
    1041             if (_delta3->state(e) == _delta3->IN_HEAP) {
    1042               _delta3->erase(e);
    1043             }
    10441032          } else {
    10451033
     
    10781066          std::numeric_limits<Value>::max()) {
    10791067        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
    1080                        (*_blossom_data)[blossom].offset);
     1068                      (*_blossom_data)[blossom].offset);
    10811069      }
    10821070
     
    11171105            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
    11181106              _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);
    11231107            }
    11241108          } else {
     
    11581142    }
    11591143
    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 
    12551144    void alternatePath(int even, int tree) {
    12561145      int odd;
     
    12951184      destroyTree(tree);
    12961185
    1297       (*_blossom_data)[blossom].status = UNMATCHED;
    12981186      (*_blossom_data)[blossom].base = node;
    1299       matchedToUnmatched(blossom);
    1300     }
    1301 
     1187      (*_blossom_data)[blossom].next = INVALID;
     1188    }
    13021189
    13031190    void augmentOnEdge(const Edge& edge) {
     
    13061193      int right = _blossom_set->find(_graph.v(edge));
    13071194
    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       }
     1195      int left_tree = _tree_set->find(left);
     1196      alternatePath(left, left_tree);
     1197      destroyTree(left_tree);
     1198
     1199      int right_tree = _tree_set->find(right);
     1200      alternatePath(right, right_tree);
     1201      destroyTree(right_tree);
    13251202
    13261203      (*_blossom_data)[left].next = _graph.direct(edge, true);
    13271204      (*_blossom_data)[right].next = _graph.direct(edge, false);
     1205    }
     1206
     1207    void augmentOnArc(const Arc& arc) {
     1208
     1209      int left = _blossom_set->find(_graph.source(arc));
     1210      int right = _blossom_set->find(_graph.target(arc));
     1211
     1212      (*_blossom_data)[left].status = MATCHED;
     1213
     1214      int right_tree = _tree_set->find(right);
     1215      alternatePath(right, right_tree);
     1216      destroyTree(right_tree);
     1217
     1218      (*_blossom_data)[left].next = arc;
     1219      (*_blossom_data)[right].next = _graph.oppositeArc(arc);
    13281220    }
    13291221
     
    15301422          (*_blossom_data)[sb].pred = pred;
    15311423          (*_blossom_data)[sb].next =
    1532                            _graph.oppositeArc((*_blossom_data)[tb].next);
     1424            _graph.oppositeArc((*_blossom_data)[tb].next);
    15331425
    15341426          pred = (*_blossom_data)[ub].next;
     
    16301522
    16311523      for (int i = 0; i < int(blossoms.size()); ++i) {
    1632         if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
     1524        if ((*_blossom_data)[blossoms[i]].next != INVALID) {
    16331525
    16341526          Value offset = (*_blossom_data)[blossoms[i]].offset;
     
    17581650          _delta4->prio() : std::numeric_limits<Value>::max();
    17591651
    1760         _delta_sum = d1; OpType ot = D1;
     1652        _delta_sum = d3; OpType ot = D3;
     1653        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
    17611654        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
    1762         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
    17631655        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
    1764 
    17651656
    17661657        switch (ot) {
     
    17761667            int blossom = _delta2->top();
    17771668            Node n = _blossom_set->classTop(blossom);
    1778             Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
    1779             extendOnArc(e);
     1669            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
     1670            if ((*_blossom_data)[blossom].next == INVALID) {
     1671              augmentOnArc(a);
     1672              --unmatched;
     1673            } else {
     1674              extendOnArc(a);
     1675            }
    17801676          }
    17811677          break;
     
    17901686              _delta3->pop();
    17911687            } 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               }
     1688              int left_tree = _tree_set->find(left_blossom);
     1689              int right_tree = _tree_set->find(right_blossom);
    18061690
    18071691              if (left_tree == right_tree) {
     
    18381722
    18391723    /// \name Primal Solution
    1840     /// Functions to get the primal solution, i.e. the maximum weighted 
     1724    /// Functions to get the primal solution, i.e. the maximum weighted
    18411725    /// matching.\n
    18421726    /// Either \ref run() or \ref start() function should be called before
     
    18571741        }
    18581742      }
    1859       return sum /= 2;
     1743      return sum / 2;
    18601744    }
    18611745
     
    18771761    /// \brief Return \c true if the given edge is in the matching.
    18781762    ///
    1879     /// This function returns \c true if the given edge is in the found 
     1763    /// This function returns \c true if the given edge is in the found
    18801764    /// matching.
    18811765    ///
     
    18881772    ///
    18891773    /// 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 
     1774    /// given node in the found matching or \c INVALID if the node is
    18911775    /// not covered by the matching.
    18921776    ///
     
    19061790    /// \brief Return the mate of the given node.
    19071791    ///
    1908     /// This function returns the mate of the given node in the found 
     1792    /// This function returns the mate of the given node in the found
    19091793    /// matching or \c INVALID if the node is not covered by the matching.
    19101794    ///
     
    19261810    /// \brief Return the value of the dual solution.
    19271811    ///
    1928     /// This function returns the value of the dual solution. 
    1929     /// It should be equal to the primal value scaled by \ref dualScale 
     1812    /// This function returns the value of the dual solution.
     1813    /// It should be equal to the primal value scaled by \ref dualScale
    19301814    /// "dual scale".
    19311815    ///
     
    19821866    /// \brief Iterator for obtaining the nodes of a blossom.
    19831867    ///
    1984     /// This class provides an iterator for obtaining the nodes of the 
     1868    /// This class provides an iterator for obtaining the nodes of the
    19851869    /// given blossom. It lists a subset of the nodes.
    1986     /// Before using this iterator, you must allocate a 
     1870    /// Before using this iterator, you must allocate a
    19871871    /// MaxWeightedMatching class and execute it.
    19881872    class BlossomIt {
     
    19931877      /// Constructor to get the nodes of the given variable.
    19941878      ///
    1995       /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 
    1996       /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 
     1879      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
     1880      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
    19971881      /// called before initializing this iterator.
    19981882      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
     
    20471931  /// \f$O(nm\log n)\f$ time complexity.
    20481932  ///
    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 
     1933  /// The maximum weighted perfect matching problem is to find a subset of
     1934  /// the edges in an undirected graph with maximum overall weight for which
    20511935  /// each node has exactly one incident edge.
    20521936  /// It can be formulated with the following linear program.
     
    20711955      \frac{\vert B \vert - 1}{2}z_B\f] */
    20721956  ///
    2073   /// The algorithm can be executed with the run() function. 
     1957  /// The algorithm can be executed with the run() function.
    20741958  /// 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. 
     1959  /// can be obtained using the query functions and the
     1960  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
     1961  /// which is able to iterate on the nodes of a blossom.
    20781962  /// If the value type is integer, then the dual solution is multiplied
    20791963  /// by \ref MaxWeightedMatching::dualScale "4".
    20801964  ///
    20811965  /// \tparam GR The undirected graph type the algorithm runs on.
    2082   /// \tparam WM The type edge weight map. The default type is 
     1966  /// \tparam WM The type edge weight map. The default type is
    20831967  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
    20841968#ifdef DOXYGEN
     
    22342118
    22352119    void destroyStructures() {
    2236       _node_num = countNodes(_graph);
    2237       _blossom_num = _node_num * 3 / 2;
    2238 
    22392120      if (_matching) {
    22402121        delete _matching;
     
    29922873          _delta4->prio() : std::numeric_limits<Value>::max();
    29932874
    2994         _delta_sum = d2; OpType ot = D2;
    2995         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
     2875        _delta_sum = d3; OpType ot = D3;
     2876        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
    29962877        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
    29972878
     
    30562937
    30572938    /// \name Primal Solution
    3058     /// Functions to get the primal solution, i.e. the maximum weighted 
     2939    /// Functions to get the primal solution, i.e. the maximum weighted
    30592940    /// perfect matching.\n
    30602941    /// Either \ref run() or \ref start() function should be called before
     
    30752956        }
    30762957      }
    3077       return sum /= 2;
     2958      return sum / 2;
    30782959    }
    30792960
    30802961    /// \brief Return \c true if the given edge is in the matching.
    30812962    ///
    3082     /// This function returns \c true if the given edge is in the found 
     2963    /// This function returns \c true if the given edge is in the found
    30832964    /// matching.
    30842965    ///
     
    30912972    ///
    30922973    /// 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 
     2974    /// given node in the found matching or \c INVALID if the node is
    30942975    /// not covered by the matching.
    30952976    ///
     
    31092990    /// \brief Return the mate of the given node.
    31102991    ///
    3111     /// This function returns the mate of the given node in the found 
     2992    /// This function returns the mate of the given node in the found
    31122993    /// matching or \c INVALID if the node is not covered by the matching.
    31132994    ///
     
    31283009    /// \brief Return the value of the dual solution.
    31293010    ///
    3130     /// This function returns the value of the dual solution. 
    3131     /// It should be equal to the primal value scaled by \ref dualScale 
     3011    /// This function returns the value of the dual solution.
     3012    /// It should be equal to the primal value scaled by \ref dualScale
    31323013    /// "dual scale".
    31333014    ///
     
    31843065    /// \brief Iterator for obtaining the nodes of a blossom.
    31853066    ///
    3186     /// This class provides an iterator for obtaining the nodes of the 
     3067    /// This class provides an iterator for obtaining the nodes of the
    31873068    /// given blossom. It lists a subset of the nodes.
    3188     /// Before using this iterator, you must allocate a 
     3069    /// Before using this iterator, you must allocate a
    31893070    /// MaxWeightedPerfectMatching class and execute it.
    31903071    class BlossomIt {
     
    31953076      /// Constructor to get the nodes of the given variable.
    31963077      ///
    3197       /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 
    3198       /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 
     3078      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
     3079      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
    31993080      /// must be called before initializing this iterator.
    32003081      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
     
    32423123} //END OF NAMESPACE LEMON
    32433124
    3244 #endif //LEMON_MAX_MATCHING_H
     3125#endif //LEMON_MATCHING_H
Note: See TracChangeset for help on using the changeset viewer.