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 @@ -28,6 +28,7 @@ #include #include #include +#include ///\ingroup matching ///\file @@ -41,7 +42,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 +70,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 +513,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 +535,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 +557,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 +571,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 +580,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 +596,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 +606,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 +649,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 +674,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 +746,7 @@ typedef RangeMap IntIntMap; enum Status { - EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2 + EVEN = -1, MATCHED = 0, ODD = 1 }; typedef HeapUnionFind BlossomSet; @@ -797,6 +798,10 @@ BinHeap *_delta4; Value _delta_sum; + int _unmatched; + + typedef MaxWeightedFractionalMatching FractionalMatching; + FractionalMatching *_fractional; void createStructures() { _node_num = countNodes(_graph); @@ -844,9 +849,6 @@ } void destroyStructures() { - _node_num = countNodes(_graph); - _blossom_num = _node_num * 3 / 2; - if (_matching) { delete _matching; } @@ -922,10 +924,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 +947,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 +959,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 +1106,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 +1188,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 +1426,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 +1526,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; @@ -1667,10 +1564,16 @@ _delta3_index(0), _delta3(0), _delta4_index(0), _delta4(0), - _delta_sum() {} + _delta_sum(), _unmatched(0), + + _fractional(0) + {} ~MaxWeightedMatching() { destroyStructures(); + if (_fractional) { + delete _fractional; + } } /// \name Execution Control @@ -1699,6 +1602,8 @@ (*_delta4_index)[i] = _delta4->PRE_HEAP; } + _unmatched = _node_num; + int index = 0; for (NodeIt n(_graph); n != INVALID; ++n) { Value max = 0; @@ -1733,18 +1638,155 @@ } } + /// \brief Initialize the algorithm with fractional matching + /// + /// This function initializes the algorithm with a fractional + /// matching. This initialization is also called jumpstart heuristic. + void fractionalInit() { + createStructures(); + + if (_fractional == 0) { + _fractional = new FractionalMatching(_graph, _weight, false); + } + _fractional->run(); + + for (ArcIt e(_graph); e != INVALID; ++e) { + (*_node_heap_index)[e] = BinHeap::PRE_HEAP; + } + for (NodeIt n(_graph); n != INVALID; ++n) { + (*_delta1_index)[n] = _delta1->PRE_HEAP; + } + for (EdgeIt e(_graph); e != INVALID; ++e) { + (*_delta3_index)[e] = _delta3->PRE_HEAP; + } + for (int i = 0; i < _blossom_num; ++i) { + (*_delta2_index)[i] = _delta2->PRE_HEAP; + (*_delta4_index)[i] = _delta4->PRE_HEAP; + } + + _unmatched = 0; + + int index = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + Value pot = _fractional->nodeValue(n); + (*_node_index)[n] = index; + (*_node_data)[index].pot = pot; + int blossom = + _blossom_set->insert(n, std::numeric_limits::max()); + + (*_blossom_data)[blossom].status = MATCHED; + (*_blossom_data)[blossom].pred = INVALID; + (*_blossom_data)[blossom].next = _fractional->matching(n); + if (_fractional->matching(n) == INVALID) { + (*_blossom_data)[blossom].base = n; + } + (*_blossom_data)[blossom].pot = 0; + (*_blossom_data)[blossom].offset = 0; + ++index; + } + + typename Graph::template NodeMap processed(_graph, false); + for (NodeIt n(_graph); n != INVALID; ++n) { + if (processed[n]) continue; + processed[n] = true; + if (_fractional->matching(n) == INVALID) continue; + int num = 1; + Node v = _graph.target(_fractional->matching(n)); + while (n != v) { + processed[v] = true; + v = _graph.target(_fractional->matching(v)); + ++num; + } + + if (num % 2 == 1) { + std::vector subblossoms(num); + + subblossoms[--num] = _blossom_set->find(n); + _delta1->push(n, _fractional->nodeValue(n)); + v = _graph.target(_fractional->matching(n)); + while (n != v) { + subblossoms[--num] = _blossom_set->find(v); + _delta1->push(v, _fractional->nodeValue(v)); + v = _graph.target(_fractional->matching(v)); + } + + int surface = + _blossom_set->join(subblossoms.begin(), subblossoms.end()); + (*_blossom_data)[surface].status = EVEN; + (*_blossom_data)[surface].pred = INVALID; + (*_blossom_data)[surface].next = INVALID; + (*_blossom_data)[surface].pot = 0; + (*_blossom_data)[surface].offset = 0; + + _tree_set->insert(surface); + ++_unmatched; + } + } + + for (EdgeIt e(_graph); e != INVALID; ++e) { + int si = (*_node_index)[_graph.u(e)]; + int sb = _blossom_set->find(_graph.u(e)); + int ti = (*_node_index)[_graph.v(e)]; + int tb = _blossom_set->find(_graph.v(e)); + if ((*_blossom_data)[sb].status == EVEN && + (*_blossom_data)[tb].status == EVEN && sb != tb) { + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - + dualScale * _weight[e]) / 2); + } + } + + for (NodeIt n(_graph); n != INVALID; ++n) { + int nb = _blossom_set->find(n); + if ((*_blossom_data)[nb].status != MATCHED) continue; + int ni = (*_node_index)[n]; + + 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) { + + int vt = _tree_set->find(vb); + + 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, e); + (*_node_data)[ni].heap.decrease(e, rw); + it->second = e; + } + } else { + (*_node_data)[ni].heap.push(e, rw); + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); + } + } + } + + if (!(*_node_data)[ni].heap.empty()) { + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); + _delta2->push(nb, _blossom_set->classPrio(nb)); + } + } + } + /// \brief Start the algorithm /// /// This function starts the algorithm. /// - /// \pre \ref init() must be called before using this function. + /// \pre \ref init() or \ref fractionalInit() must be called + /// before using this function. void start() { enum OpType { D1, D2, D3, D4 }; - int unmatched = _node_num; - while (unmatched > 0) { + while (_unmatched > 0) { Value d1 = !_delta1->empty() ? _delta1->prio() : std::numeric_limits::max(); @@ -1757,26 +1799,30 @@ 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: { Node n = _delta1->top(); unmatchNode(n); - --unmatched; + --_unmatched; } break; case D2: { 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,26 +1835,14 @@ 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); } else { augmentOnEdge(e); - unmatched -= 2; + _unmatched -= 2; } } } break; @@ -1826,18 +1860,18 @@ /// /// \note mwm.run() is just a shortcut of the following code. /// \code - /// mwm.init(); + /// mwm.fractionalInit(); /// mwm.start(); /// \endcode void run() { - init(); + fractionalInit(); start(); } /// @} /// \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 +1890,7 @@ sum += _weight[(*_matching)[n]]; } } - return sum /= 2; + return sum / 2; } /// \brief Return the size (cardinality) of the matching. @@ -1876,7 +1910,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 +1921,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 +1939,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 +1959,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 +2015,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 +2026,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 +2080,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 +2104,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 @@ -2190,6 +2224,11 @@ BinHeap *_delta4; Value _delta_sum; + int _unmatched; + + typedef MaxWeightedPerfectFractionalMatching + FractionalMatching; + FractionalMatching *_fractional; void createStructures() { _node_num = countNodes(_graph); @@ -2233,9 +2272,6 @@ } void destroyStructures() { - _node_num = countNodes(_graph); - _blossom_num = _node_num * 3 / 2; - if (_matching) { delete _matching; } @@ -2908,10 +2944,16 @@ _delta3_index(0), _delta3(0), _delta4_index(0), _delta4(0), - _delta_sum() {} + _delta_sum(), _unmatched(0), + + _fractional(0) + {} ~MaxWeightedPerfectMatching() { destroyStructures(); + if (_fractional) { + delete _fractional; + } } /// \name Execution Control @@ -2937,6 +2979,8 @@ (*_delta4_index)[i] = _delta4->PRE_HEAP; } + _unmatched = _node_num; + int index = 0; for (NodeIt n(_graph); n != INVALID; ++n) { Value max = - std::numeric_limits::max(); @@ -2970,18 +3014,152 @@ } } + /// \brief Initialize the algorithm with fractional matching + /// + /// This function initializes the algorithm with a fractional + /// matching. This initialization is also called jumpstart heuristic. + void fractionalInit() { + createStructures(); + + if (_fractional == 0) { + _fractional = new FractionalMatching(_graph, _weight, false); + } + if (!_fractional->run()) { + _unmatched = -1; + return; + } + + for (ArcIt e(_graph); e != INVALID; ++e) { + (*_node_heap_index)[e] = BinHeap::PRE_HEAP; + } + for (EdgeIt e(_graph); e != INVALID; ++e) { + (*_delta3_index)[e] = _delta3->PRE_HEAP; + } + for (int i = 0; i < _blossom_num; ++i) { + (*_delta2_index)[i] = _delta2->PRE_HEAP; + (*_delta4_index)[i] = _delta4->PRE_HEAP; + } + + _unmatched = 0; + + int index = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + Value pot = _fractional->nodeValue(n); + (*_node_index)[n] = index; + (*_node_data)[index].pot = pot; + int blossom = + _blossom_set->insert(n, std::numeric_limits::max()); + + (*_blossom_data)[blossom].status = MATCHED; + (*_blossom_data)[blossom].pred = INVALID; + (*_blossom_data)[blossom].next = _fractional->matching(n); + (*_blossom_data)[blossom].pot = 0; + (*_blossom_data)[blossom].offset = 0; + ++index; + } + + typename Graph::template NodeMap processed(_graph, false); + for (NodeIt n(_graph); n != INVALID; ++n) { + if (processed[n]) continue; + processed[n] = true; + if (_fractional->matching(n) == INVALID) continue; + int num = 1; + Node v = _graph.target(_fractional->matching(n)); + while (n != v) { + processed[v] = true; + v = _graph.target(_fractional->matching(v)); + ++num; + } + + if (num % 2 == 1) { + std::vector subblossoms(num); + + subblossoms[--num] = _blossom_set->find(n); + v = _graph.target(_fractional->matching(n)); + while (n != v) { + subblossoms[--num] = _blossom_set->find(v); + v = _graph.target(_fractional->matching(v)); + } + + int surface = + _blossom_set->join(subblossoms.begin(), subblossoms.end()); + (*_blossom_data)[surface].status = EVEN; + (*_blossom_data)[surface].pred = INVALID; + (*_blossom_data)[surface].next = INVALID; + (*_blossom_data)[surface].pot = 0; + (*_blossom_data)[surface].offset = 0; + + _tree_set->insert(surface); + ++_unmatched; + } + } + + for (EdgeIt e(_graph); e != INVALID; ++e) { + int si = (*_node_index)[_graph.u(e)]; + int sb = _blossom_set->find(_graph.u(e)); + int ti = (*_node_index)[_graph.v(e)]; + int tb = _blossom_set->find(_graph.v(e)); + if ((*_blossom_data)[sb].status == EVEN && + (*_blossom_data)[tb].status == EVEN && sb != tb) { + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - + dualScale * _weight[e]) / 2); + } + } + + for (NodeIt n(_graph); n != INVALID; ++n) { + int nb = _blossom_set->find(n); + if ((*_blossom_data)[nb].status != MATCHED) continue; + int ni = (*_node_index)[n]; + + 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) { + + int vt = _tree_set->find(vb); + + 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, e); + (*_node_data)[ni].heap.decrease(e, rw); + it->second = e; + } + } else { + (*_node_data)[ni].heap.push(e, rw); + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); + } + } + } + + if (!(*_node_data)[ni].heap.empty()) { + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); + _delta2->push(nb, _blossom_set->classPrio(nb)); + } + } + } + /// \brief Start the algorithm /// /// This function starts the algorithm. /// - /// \pre \ref init() must be called before using this function. + /// \pre \ref init() or \ref fractionalInit() must be called before + /// using this function. bool start() { enum OpType { D2, D3, D4 }; - int unmatched = _node_num; - while (unmatched > 0) { + if (_unmatched == -1) return false; + + while (_unmatched > 0) { Value d2 = !_delta2->empty() ? _delta2->prio() : std::numeric_limits::max(); @@ -2991,8 +3169,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()) { @@ -3025,7 +3203,7 @@ shrinkOnEdge(e, left_tree); } else { augmentOnEdge(e); - unmatched -= 2; + _unmatched -= 2; } } } break; @@ -3044,18 +3222,18 @@ /// /// \note mwpm.run() is just a shortcut of the following code. /// \code - /// mwpm.init(); + /// mwpm.fractionalInit(); /// mwpm.start(); /// \endcode bool run() { - init(); + fractionalInit(); return start(); } /// @} /// \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 +3252,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 +3268,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 +3286,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 +3305,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 +3361,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 +3372,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 +3419,4 @@ } //END OF NAMESPACE LEMON -#endif //LEMON_MAX_MATCHING_H +#endif //LEMON_MATCHING_H