# HG changeset patch # User Balazs Dezso # Date 2009-09-20 21:38:24 # Node ID 0513ccfea9670a94312cc2bde2ed8595f08361aa # Parent 6d5f547e5bfb8131568ad50ea406129a056ac9ed 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. diff --git a/lemon/matching.h b/lemon/matching.h --- a/lemon/matching.h +++ b/lemon/matching.h @@ -16,8 +16,8 @@ * */ -#ifndef LEMON_MAX_MATCHING_H -#define LEMON_MAX_MATCHING_H +#ifndef LEMON_MATCHING_H +#define LEMON_MATCHING_H #include #include @@ -41,7 +41,7 @@ /// /// 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 @@ -69,11 +69,11 @@ ///\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.) @@ -512,7 +512,7 @@ } } - /// \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 @@ -534,8 +534,8 @@ /// \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)) { @@ -556,7 +556,7 @@ /// \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; @@ -570,7 +570,7 @@ /// \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)]; @@ -579,7 +579,7 @@ /// \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]; @@ -595,7 +595,7 @@ /// \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 ? @@ -605,7 +605,7 @@ /// @} /// \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. /// @{ @@ -648,8 +648,8 @@ /// 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] @@ -673,16 +673,16 @@ /** \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 @@ -745,7 +745,7 @@ typedef RangeMap IntIntMap; enum Status { - EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2 + EVEN = -1, MATCHED = 0, ODD = 1 }; typedef HeapUnionFind BlossomSet; @@ -844,9 +844,6 @@ } void destroyStructures() { - _node_num = countNodes(_graph); - _blossom_num = _node_num * 3 / 2; - if (_matching) { delete _matching; } @@ -922,10 +919,6 @@ 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); @@ -949,202 +942,6 @@ _delta2->push(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset){ - _delta2->decrease(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } - } - } - } - } - } - (*_blossom_data)[blossom].offset = 0; - } - - void matchedToOdd(int blossom) { - if (_delta2->state(blossom) == _delta2->IN_HEAP) { - _delta2->erase(blossom); - } - (*_blossom_data)[blossom].offset += _delta_sum; - if (!_blossom_set->trivial(blossom)) { - _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + - (*_blossom_data)[blossom].offset); - } - } - - void evenToMatched(int blossom, int tree) { - if (!_blossom_set->trivial(blossom)) { - (*_blossom_data)[blossom].pot += 2 * _delta_sum; - } - - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); - n != INVALID; ++n) { - int ni = (*_node_index)[n]; - (*_node_data)[ni].pot -= _delta_sum; - - _delta1->erase(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); - - if (vt != tree) { - - 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); - } - } else { - - typename std::map::iterator it = - (*_node_data)[vi].heap_index.find(tree); - - if (it != (*_node_data)[vi].heap_index.end()) { - (*_node_data)[vi].heap.erase(it->second); - (*_node_data)[vi].heap_index.erase(it); - if ((*_node_data)[vi].heap.empty()) { - _blossom_set->increase(v, std::numeric_limits::max()); - } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { - _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); - } - - if ((*_blossom_data)[vb].status == MATCHED) { - if (_blossom_set->classPrio(vb) == - std::numeric_limits::max()) { - _delta2->erase(vb); - } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset) { - _delta2->increase(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } - } - } - } - } - } - } - - void oddToMatched(int blossom) { - (*_blossom_data)[blossom].offset -= _delta_sum; - - if (_blossom_set->classPrio(blossom) != - std::numeric_limits::max()) { - _delta2->push(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); - } - - if (!_blossom_set->trivial(blossom)) { - _delta4->erase(blossom); - } - } - - void oddToEven(int blossom, int tree) { - if (!_blossom_set->trivial(blossom)) { - _delta4->erase(blossom); - (*_blossom_data)[blossom].pot -= - 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); - } - - 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(); - (*_node_data)[ni].pot += - 2 * _delta_sum - (*_blossom_data)[blossom].offset; - - _delta1->push(n, (*_node_data)[ni].pot); - - 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 ((*_blossom_data)[vb].status == EVEN) { - 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); - - if (it != (*_node_data)[vi].heap_index.end()) { - if ((*_node_data)[vi].heap[it->second] > rw) { - (*_node_data)[vi].heap.replace(it->second, e); - (*_node_data)[vi].heap.decrease(e, rw); - it->second = e; - } - } else { - (*_node_data)[vi].heap.push(e, rw); - (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); - } - - if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { - _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); - - if ((*_blossom_data)[vb].status == MATCHED) { - if (_delta2->state(vb) != _delta2->IN_HEAP) { - _delta2->push(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset) { _delta2->decrease(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); @@ -1157,43 +954,145 @@ (*_blossom_data)[blossom].offset = 0; } - - void matchedToUnmatched(int blossom) { + void matchedToOdd(int blossom) { if (_delta2->state(blossom) == _delta2->IN_HEAP) { _delta2->erase(blossom); } + (*_blossom_data)[blossom].offset += _delta_sum; + if (!_blossom_set->trivial(blossom)) { + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + + (*_blossom_data)[blossom].offset); + } + } + + void evenToMatched(int blossom, int tree) { + if (!_blossom_set->trivial(blossom)) { + (*_blossom_data)[blossom].pot += 2 * _delta_sum; + } 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); + (*_node_data)[ni].pot -= _delta_sum; + + _delta1->erase(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 ((*_blossom_data)[vb].status == EVEN) { - if (_delta3->state(e) != _delta3->IN_HEAP) { - _delta3->push(e, rw); + 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); + + if (vt != tree) { + + 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 { + + typename std::map::iterator it = + (*_node_data)[vi].heap_index.find(tree); + + if (it != (*_node_data)[vi].heap_index.end()) { + (*_node_data)[vi].heap.erase(it->second); + (*_node_data)[vi].heap_index.erase(it); + if ((*_node_data)[vi].heap.empty()) { + _blossom_set->increase(v, std::numeric_limits::max()); + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { + _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); + } + + if ((*_blossom_data)[vb].status == MATCHED) { + if (_blossom_set->classPrio(vb) == + std::numeric_limits::max()) { + _delta2->erase(vb); + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset) { + _delta2->increase(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } + } } } } } } - void unmatchedToMatched(int blossom) { + void oddToMatched(int blossom) { + (*_blossom_data)[blossom].offset -= _delta_sum; + + if (_blossom_set->classPrio(blossom) != + std::numeric_limits::max()) { + _delta2->push(blossom, _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset); + } + + if (!_blossom_set->trivial(blossom)) { + _delta4->erase(blossom); + } + } + + void oddToEven(int blossom, int tree) { + if (!_blossom_set->trivial(blossom)) { + _delta4->erase(blossom); + (*_blossom_data)[blossom].pot -= + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); + } + 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(); + (*_node_data)[ni].pot += + 2 * _delta_sum - (*_blossom_data)[blossom].offset; + + _delta1->push(n, (*_node_data)[ni].pot); + for (InArcIt e(_graph, n); e != INVALID; ++e) { Node v = _graph.source(e); int vb = _blossom_set->find(v); @@ -1202,54 +1101,44 @@ 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); + if ((*_blossom_data)[vb].status == EVEN) { + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { + _delta3->push(e, rw / 2); } - } 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); + } else { 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; + (*_node_data)[vi].heap_index.find(tree); + + if (it != (*_node_data)[vi].heap_index.end()) { + if ((*_node_data)[vi].heap[it->second] > rw) { + (*_node_data)[vi].heap.replace(it->second, e); + (*_node_data)[vi].heap.decrease(e, rw); + it->second = e; } } else { - (*_node_data)[ni].heap.push(r, rw); - (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); + (*_node_data)[vi].heap.push(e, rw); + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); } - 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); + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); + + if ((*_blossom_data)[vb].status == MATCHED) { + if (_delta2->state(vb) != _delta2->IN_HEAP) { + _delta2->push(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset) { + _delta2->decrease(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } } } - - } else if ((*_blossom_data)[vb].status == UNMATCHED) { - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } } } } + (*_blossom_data)[blossom].offset = 0; } void alternatePath(int even, int tree) { @@ -1294,39 +1183,42 @@ 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); - } - - 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); + } + void extendOnArc(const Arc& arc) { int base = _blossom_set->find(_graph.target(arc)); int tree = _tree_set->find(base); @@ -1529,7 +1421,7 @@ _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; @@ -1629,7 +1521,7 @@ } 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; @@ -1757,12 +1649,11 @@ 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: { @@ -1775,8 +1666,13 @@ { 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: @@ -1789,20 +1685,8 @@ 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); @@ -1837,7 +1721,7 @@ /// @} /// \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. @@ -1856,7 +1740,7 @@ sum += _weight[(*_matching)[n]]; } } - return sum /= 2; + return sum / 2; } /// \brief Return the size (cardinality) of the matching. @@ -1876,7 +1760,7 @@ /// \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. @@ -1887,7 +1771,7 @@ /// \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. @@ -1905,7 +1789,7 @@ /// \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. @@ -1925,8 +1809,8 @@ /// \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. @@ -1981,9 +1865,9 @@ /// \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: @@ -1992,8 +1876,8 @@ /// /// 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) @@ -2046,8 +1930,8 @@ /// 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] @@ -2070,16 +1954,16 @@ /** \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 @@ -2233,9 +2117,6 @@ } void destroyStructures() { - _node_num = countNodes(_graph); - _blossom_num = _node_num * 3 / 2; - if (_matching) { delete _matching; } @@ -2991,8 +2872,8 @@ 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()) { @@ -3055,7 +2936,7 @@ /// @} /// \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. @@ -3074,12 +2955,12 @@ 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. @@ -3090,7 +2971,7 @@ /// \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. @@ -3108,7 +2989,7 @@ /// \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. @@ -3127,8 +3008,8 @@ /// \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. @@ -3183,9 +3064,9 @@ /// \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: @@ -3194,8 +3075,8 @@ /// /// 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) @@ -3241,4 +3122,4 @@ } //END OF NAMESPACE LEMON -#endif //LEMON_MAX_MATCHING_H +#endif //LEMON_MATCHING_H