# Ticket #314: 0513ccfea967.patch

File 0513ccfea967.patch, 25.5 KB (added by Balazs Dezso, 15 years ago)

Improvements in matching.h

• ## lemon/matching.h

# HG changeset patch
# User Balazs Dezso <deba@inf.elte.hu>
# Date 1253475504 -7200
# Node ID 0513ccfea9670a94312cc2bde2ed8595f08361aa
diff -r 6d5f547e5bfb -r 0513ccfea967 lemon/matching.h
 a * */ #ifndef LEMON_MAX_MATCHING_H #define LEMON_MAX_MATCHING_H #ifndef LEMON_MATCHING_H #define LEMON_MATCHING_H #include #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). /// /// The dual solution of the problem is a map of the nodes to ///\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 { EVEN = 1,       ///< = 1. (\c D is an alias for \c EVEN.) } } /// \brief Start Edmonds' algorithm with a heuristic improvement /// \brief Start Edmonds' algorithm with a heuristic improvement /// for dense graphs /// /// This function runs Edmonds' algorithm with a heuristic of postponing /// \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() { if (countEdges(_graph) < 2 * countNodes(_graph)) { /// \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 { int size = 0; /// \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 { return edge == (*_matching)[_graph.u(edge)]; /// \brief Return the matching arc (or edge) incident to the given node. /// /// 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 { return (*_matching)[n]; /// \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 { return (*_matching)[n] != INVALID ? /// @} /// \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. /// @{ /// on extensive use of priority queues and provides /// \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. /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f] /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} \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 template typedef RangeMap IntIntMap; enum Status { EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2 EVEN = -1, MATCHED = 0, ODD = 1 }; typedef HeapUnionFind BlossomSet; } 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 { typename std::map::iterator it = (*_node_data)[vi].heap_index.find(tree); _delta2->push(vb, _blossom_set->classPrio(vb) - (*_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); } (*_blossom_data)[blossom].offset += _delta_sum; 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 { typename std::map::iterator it = if (_blossom_set->classPrio(blossom) != std::numeric_limits::max()) { _delta2->push(blossom, _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset); (*_blossom_data)[blossom].offset); } if (!_blossom_set->trivial(blossom)) { 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 { typename std::map::iterator it = (*_blossom_data)[blossom].offset = 0; } 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; alternatePath(blossom, tree); 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 left = _blossom_set->find(_graph.u(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); } int left_tree = _tree_set->find(left); alternatePath(left, left_tree); destroyTree(left_tree); 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 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); } void extendOnArc(const Arc& arc) { int base = _blossom_set->find(_graph.target(arc)); int tree = _tree_set->find(base); _tree_set->insert(sb, tree); (*_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; (*_blossom_data)[blossoms[i]].pot += 2 * offset; Value d4 = !_delta4->empty() ? _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) { case D1: { { 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; case D3: if (left_blossom == right_blossom) { _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) { shrinkOnEdge(e, left_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 /// using them. sum += _weight[(*_matching)[n]]; } } return sum /= 2; return sum / 2; } /// \brief Return the size (cardinality) of the matching. /// \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. /// /// \pre Either run() or start() must be called before using this function. /// \brief Return the matching arc (or edge) incident to the given node. /// /// 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. /// /// \pre Either run() or start() must be called before using this function. /// \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. /// /// \pre Either run() or start() must be called before using this function. /// \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". /// /// \pre Either run() or start() must be called before using this function. /// \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 { public: /// /// 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) : _algorithm(&algorithm) /// is based on extensive use of priority queues and provides /// \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. /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f] /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} \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 template } void destroyStructures() { _node_num = countNodes(_graph); _blossom_num = _node_num * 3 / 2; if (_matching) { delete _matching; } Value d4 = !_delta4->empty() ? _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; } if (_delta_sum == std::numeric_limits::max()) { /// @} /// \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 /// using them. sum += _weight[(*_matching)[n]]; } } 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. /// /// \pre Either run() or start() must be called before using this function. /// \brief Return the matching arc (or edge) incident to the given node. /// /// 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. /// /// \pre Either run() or start() must be called before using this function. /// \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. /// /// \pre Either run() or start() must be called before using this function. /// \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". /// /// \pre Either run() or start() must be called before using this function. /// \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 { public: /// /// 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) : _algorithm(&algorithm) } //END OF NAMESPACE LEMON #endif //LEMON_MAX_MATCHING_H #endif //LEMON_MATCHING_H