# Changeset 868:0513ccfea967 in lemon-1.2 for lemon

Ignore:
Timestamp:
09/20/09 21:38:24 (10 years ago)
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
 r651 */ #ifndef LEMON_MAX_MATCHING_H #define LEMON_MAX_MATCHING_H #ifndef LEMON_MATCHING_H #define LEMON_MATCHING_H #include /// This class implements Edmonds' alternating forest matching algorithm /// for finding a maximum cardinality matching in a general undirected graph. /// It can be started from an arbitrary initial matching /// It can be started from an arbitrary initial matching /// (the default is the empty one). /// ///\brief Status constants for Gallai-Edmonds decomposition. /// ///These constants are used for indicating the Gallai-Edmonds ///These constants are used for indicating the Gallai-Edmonds ///decomposition of a graph. The nodes with status \c EVEN (or \c D) ///induce a subgraph with factor-critical components, the nodes with ///status \c ODD (or \c A) form the canonical barrier, and the nodes ///with status \c MATCHED (or \c C) induce a subgraph having a ///with status \c MATCHED (or \c C) induce a subgraph having a ///perfect matching. enum Status { } /// \brief Start Edmonds' algorithm with a heuristic improvement /// \brief Start Edmonds' algorithm with a heuristic improvement /// for dense graphs /// /// \brief Run Edmonds' algorithm /// /// This function runs Edmonds' algorithm. An additional heuristic of /// postponing shrinks is used for relatively dense graphs /// This function runs Edmonds' algorithm. An additional heuristic of /// postponing shrinks is used for relatively dense graphs /// (for which m>=2*n holds). void run() { /// \brief Return the size (cardinality) of the matching. /// /// This function returns the size (cardinality) of the current matching. /// This function returns the size (cardinality) of the current matching. /// After run() it returns the size of the maximum matching in the graph. int matchingSize() const { /// \brief Return \c true if the given edge is in the matching. /// /// This function returns \c true if the given edge is in the current /// This function returns \c true if the given edge is in the current /// matching. bool matching(const Edge& edge) const { /// /// This function returns the matching arc (or edge) incident to the /// given node in the current matching or \c INVALID if the node is /// given node in the current matching or \c INVALID if the node is /// not covered by the matching. Arc matching(const Node& n) const { /// \brief Return the mate of the given node. /// /// This function returns the mate of the given node in the current /// This function returns the mate of the given node in the current /// matching or \c INVALID if the node is not covered by the matching. Node mate(const Node& n) const { /// \name Dual Solution /// Functions to get the dual solution, i.e. the Gallai-Edmonds /// Functions to get the dual solution, i.e. the Gallai-Edmonds /// decomposition. /// \f$O(nm\log n)\f$ time complexity. /// /// The maximum weighted matching problem is to find a subset of the /// edges in an undirected graph with maximum overall weight for which /// The maximum weighted matching problem is to find a subset of the /// edges in an undirected graph with maximum overall weight for which /// each node has at most one incident edge. /// It can be formulated with the following linear program. \frac{\vert B \vert - 1}{2}z_B\f] */ /// /// The algorithm can be executed with the run() function. /// The algorithm can be executed with the run() function. /// After it the matching (the primal solution) and the dual solution /// can be obtained using the query functions and the /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, /// which is able to iterate on the nodes of a blossom. /// can be obtained using the query functions and the /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, /// which is able to iterate on the nodes of a blossom. /// If the value type is integer, then the dual solution is multiplied /// by \ref MaxWeightedMatching::dualScale "4". /// /// \tparam GR The undirected graph type the algorithm runs on. /// \tparam WM The type edge weight map. The default type is /// \tparam WM The type edge weight map. The default type is /// \ref concepts::Graph::EdgeMap "GR::EdgeMap". #ifdef DOXYGEN enum Status { EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2 EVEN = -1, MATCHED = 0, ODD = 1 }; void destroyStructures() { _node_num = countNodes(_graph); _blossom_num = _node_num * 3 / 2; if (_matching) { delete _matching; if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { _delta3->push(e, rw / 2); } } else if ((*_blossom_data)[vb].status == UNMATCHED) { if (_delta3->state(e) != _delta3->IN_HEAP) { _delta3->push(e, rw); } } else { (*_blossom_data)[vb].offset); } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset){ (*_blossom_data)[vb].offset) { _delta2->decrease(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); if (!_blossom_set->trivial(blossom)) { _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + (*_blossom_data)[blossom].offset); (*_blossom_data)[blossom].offset); } } } } } else if ((*_blossom_data)[vb].status == UNMATCHED) { if (_delta3->state(e) == _delta3->IN_HEAP) { _delta3->erase(e); } } else { std::numeric_limits::max()) { _delta2->push(blossom, _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset); (*_blossom_data)[blossom].offset); } if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { _delta3->push(e, rw / 2); } } else if ((*_blossom_data)[vb].status == UNMATCHED) { if (_delta3->state(e) != _delta3->IN_HEAP) { _delta3->push(e, rw); } } else { } void matchedToUnmatched(int blossom) { if (_delta2->state(blossom) == _delta2->IN_HEAP) { _delta2->erase(blossom); } for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); n != INVALID; ++n) { int ni = (*_node_index)[n]; _blossom_set->increase(n, std::numeric_limits::max()); (*_node_data)[ni].heap.clear(); (*_node_data)[ni].heap_index.clear(); for (OutArcIt e(_graph, n); e != INVALID; ++e) { Node v = _graph.target(e); int vb = _blossom_set->find(v); int vi = (*_node_index)[v]; Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - dualScale * _weight[e]; if ((*_blossom_data)[vb].status == EVEN) { if (_delta3->state(e) != _delta3->IN_HEAP) { _delta3->push(e, rw); } } } } } void unmatchedToMatched(int blossom) { for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); n != INVALID; ++n) { int ni = (*_node_index)[n]; for (InArcIt e(_graph, n); e != INVALID; ++e) { Node v = _graph.source(e); int vb = _blossom_set->find(v); int vi = (*_node_index)[v]; Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - dualScale * _weight[e]; if (vb == blossom) { if (_delta3->state(e) == _delta3->IN_HEAP) { _delta3->erase(e); } } else if ((*_blossom_data)[vb].status == EVEN) { if (_delta3->state(e) == _delta3->IN_HEAP) { _delta3->erase(e); } int vt = _tree_set->find(vb); Arc r = _graph.oppositeArc(e); typename std::map::iterator it = (*_node_data)[ni].heap_index.find(vt); if (it != (*_node_data)[ni].heap_index.end()) { if ((*_node_data)[ni].heap[it->second] > rw) { (*_node_data)[ni].heap.replace(it->second, r); (*_node_data)[ni].heap.decrease(r, rw); it->second = r; } } else { (*_node_data)[ni].heap.push(r, rw); (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); } if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); if (_delta2->state(blossom) != _delta2->IN_HEAP) { _delta2->push(blossom, _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset); } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)- (*_blossom_data)[blossom].offset){ _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset); } } } else if ((*_blossom_data)[vb].status == UNMATCHED) { if (_delta3->state(e) == _delta3->IN_HEAP) { _delta3->erase(e); } } } } } void alternatePath(int even, int tree) { int odd; destroyTree(tree); (*_blossom_data)[blossom].status = UNMATCHED; (*_blossom_data)[blossom].base = node; matchedToUnmatched(blossom); } (*_blossom_data)[blossom].next = INVALID; } void augmentOnEdge(const Edge& edge) { int right = _blossom_set->find(_graph.v(edge)); if ((*_blossom_data)[left].status == EVEN) { int left_tree = _tree_set->find(left); alternatePath(left, left_tree); destroyTree(left_tree); } else { (*_blossom_data)[left].status = MATCHED; unmatchedToMatched(left); } if ((*_blossom_data)[right].status == EVEN) { int right_tree = _tree_set->find(right); alternatePath(right, right_tree); destroyTree(right_tree); } else { (*_blossom_data)[right].status = MATCHED; unmatchedToMatched(right); } int left_tree = _tree_set->find(left); alternatePath(left, left_tree); destroyTree(left_tree); int right_tree = _tree_set->find(right); alternatePath(right, right_tree); destroyTree(right_tree); (*_blossom_data)[left].next = _graph.direct(edge, true); (*_blossom_data)[right].next = _graph.direct(edge, false); } void augmentOnArc(const Arc& arc) { int left = _blossom_set->find(_graph.source(arc)); int right = _blossom_set->find(_graph.target(arc)); (*_blossom_data)[left].status = MATCHED; int right_tree = _tree_set->find(right); alternatePath(right, right_tree); destroyTree(right_tree); (*_blossom_data)[left].next = arc; (*_blossom_data)[right].next = _graph.oppositeArc(arc); } (*_blossom_data)[sb].pred = pred; (*_blossom_data)[sb].next = _graph.oppositeArc((*_blossom_data)[tb].next); _graph.oppositeArc((*_blossom_data)[tb].next); pred = (*_blossom_data)[ub].next; for (int i = 0; i < int(blossoms.size()); ++i) { if ((*_blossom_data)[blossoms[i]].status == MATCHED) { if ((*_blossom_data)[blossoms[i]].next != INVALID) { Value offset = (*_blossom_data)[blossoms[i]].offset; _delta4->prio() : std::numeric_limits::max(); _delta_sum = d1; OpType ot = D1; _delta_sum = d3; OpType ot = D3; if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; } if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } switch (ot) { int blossom = _delta2->top(); Node n = _blossom_set->classTop(blossom); Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); extendOnArc(e); Arc a = (*_node_data)[(*_node_index)[n]].heap.top(); if ((*_blossom_data)[blossom].next == INVALID) { augmentOnArc(a); --unmatched; } else { extendOnArc(a); } } break; _delta3->pop(); } else { int left_tree; if ((*_blossom_data)[left_blossom].status == EVEN) { left_tree = _tree_set->find(left_blossom); } else { left_tree = -1; ++unmatched; } int right_tree; if ((*_blossom_data)[right_blossom].status == EVEN) { right_tree = _tree_set->find(right_blossom); } else { right_tree = -1; ++unmatched; } int left_tree = _tree_set->find(left_blossom); int right_tree = _tree_set->find(right_blossom); if (left_tree == right_tree) { /// \name Primal Solution /// Functions to get the primal solution, i.e. the maximum weighted /// Functions to get the primal solution, i.e. the maximum weighted /// matching.\n /// Either \ref run() or \ref start() function should be called before } } return sum /= 2; return sum / 2; } /// \brief Return \c true if the given edge is in the matching. /// /// This function returns \c true if the given edge is in the found /// This function returns \c true if the given edge is in the found /// matching. /// /// /// This function returns the matching arc (or edge) incident to the /// given node in the found matching or \c INVALID if the node is /// given node in the found matching or \c INVALID if the node is /// not covered by the matching. /// /// \brief Return the mate of the given node. /// /// This function returns the mate of the given node in the found /// This function returns the mate of the given node in the found /// matching or \c INVALID if the node is not covered by the matching. /// /// \brief Return the value of the dual solution. /// /// This function returns the value of the dual solution. /// It should be equal to the primal value scaled by \ref dualScale /// This function returns the value of the dual solution. /// It should be equal to the primal value scaled by \ref dualScale /// "dual scale". /// /// \brief Iterator for obtaining the nodes of a blossom. /// /// This class provides an iterator for obtaining the nodes of the /// This class provides an iterator for obtaining the nodes of the /// given blossom. It lists a subset of the nodes. /// Before using this iterator, you must allocate a /// Before using this iterator, you must allocate a /// MaxWeightedMatching class and execute it. class BlossomIt { /// Constructor to get the nodes of the given variable. /// /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or /// \ref MaxWeightedMatching::start() "algorithm.start()" must be /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or /// \ref MaxWeightedMatching::start() "algorithm.start()" must be /// called before initializing this iterator. BlossomIt(const MaxWeightedMatching& algorithm, int variable) /// \f$O(nm\log n)\f$ time complexity. /// /// The maximum weighted perfect matching problem is to find a subset of /// the edges in an undirected graph with maximum overall weight for which /// The maximum weighted perfect matching problem is to find a subset of /// the edges in an undirected graph with maximum overall weight for which /// each node has exactly one incident edge. /// It can be formulated with the following linear program. \frac{\vert B \vert - 1}{2}z_B\f] */ /// /// The algorithm can be executed with the run() function. /// The algorithm can be executed with the run() function. /// After it the matching (the primal solution) and the dual solution /// can be obtained using the query functions and the /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, /// which is able to iterate on the nodes of a blossom. /// can be obtained using the query functions and the /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, /// which is able to iterate on the nodes of a blossom. /// If the value type is integer, then the dual solution is multiplied /// by \ref MaxWeightedMatching::dualScale "4". /// /// \tparam GR The undirected graph type the algorithm runs on. /// \tparam WM The type edge weight map. The default type is /// \tparam WM The type edge weight map. The default type is /// \ref concepts::Graph::EdgeMap "GR::EdgeMap". #ifdef DOXYGEN void destroyStructures() { _node_num = countNodes(_graph); _blossom_num = _node_num * 3 / 2; if (_matching) { delete _matching; _delta4->prio() : std::numeric_limits::max(); _delta_sum = d2; OpType ot = D2; if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } _delta_sum = d3; OpType ot = D3; if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } /// \name Primal Solution /// Functions to get the primal solution, i.e. the maximum weighted /// Functions to get the primal solution, i.e. the maximum weighted /// perfect matching.\n /// Either \ref run() or \ref start() function should be called before } } return sum /= 2; return sum / 2; } /// \brief Return \c true if the given edge is in the matching. /// /// This function returns \c true if the given edge is in the found /// This function returns \c true if the given edge is in the found /// matching. /// /// /// This function returns the matching arc (or edge) incident to the /// given node in the found matching or \c INVALID if the node is /// given node in the found matching or \c INVALID if the node is /// not covered by the matching. /// /// \brief Return the mate of the given node. /// /// This function returns the mate of the given node in the found /// This function returns the mate of the given node in the found /// matching or \c INVALID if the node is not covered by the matching. /// /// \brief Return the value of the dual solution. /// /// This function returns the value of the dual solution. /// It should be equal to the primal value scaled by \ref dualScale /// This function returns the value of the dual solution. /// It should be equal to the primal value scaled by \ref dualScale /// "dual scale". /// /// \brief Iterator for obtaining the nodes of a blossom. /// /// This class provides an iterator for obtaining the nodes of the /// This class provides an iterator for obtaining the nodes of the /// given blossom. It lists a subset of the nodes. /// Before using this iterator, you must allocate a /// Before using this iterator, you must allocate a /// MaxWeightedPerfectMatching class and execute it. class BlossomIt { /// Constructor to get the nodes of the given variable. /// /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" /// must be called before initializing this iterator. BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) } //END OF NAMESPACE LEMON #endif //LEMON_MAX_MATCHING_H #endif //LEMON_MATCHING_H