deba@326: /* -*- mode: C++; indent-tabs-mode: nil; -*- deba@326: * deba@326: * This file is a part of LEMON, a generic C++ optimization library. deba@326: * alpar@1092: * Copyright (C) 2003-2013 deba@326: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport deba@326: * (Egervary Research Group on Combinatorial Optimization, EGRES). deba@326: * deba@326: * Permission to use, modify and distribute this software is granted deba@326: * provided that this copyright notice appears in all copies. For deba@326: * precise terms see the accompanying LICENSE file. deba@326: * deba@326: * This software is provided "AS IS" with no warranty of any kind, deba@326: * express or implied, and with no claim as to its suitability for any deba@326: * purpose. deba@326: * deba@326: */ deba@326: deba@868: #ifndef LEMON_MATCHING_H deba@868: #define LEMON_MATCHING_H deba@326: deba@326: #include deba@326: #include deba@326: #include deba@326: #include deba@326: deba@326: #include deba@326: #include deba@326: #include deba@326: #include deba@870: #include deba@326: deba@326: ///\ingroup matching deba@326: ///\file deba@327: ///\brief Maximum matching algorithms in general graphs. deba@326: deba@326: namespace lemon { deba@326: deba@327: /// \ingroup matching deba@326: /// kpeter@590: /// \brief Maximum cardinality matching in general graphs deba@326: /// kpeter@590: /// This class implements Edmonds' alternating forest matching algorithm kpeter@593: /// for finding a maximum cardinality matching in a general undirected graph. deba@868: /// It can be started from an arbitrary initial matching kpeter@590: /// (the default is the empty one). deba@326: /// alpar@330: /// The dual solution of the problem is a map of the nodes to kpeter@590: /// \ref MaxMatching::Status "Status", having values \c EVEN (or \c D), kpeter@590: /// \c ODD (or \c A) and \c MATCHED (or \c C) defining the Gallai-Edmonds kpeter@590: /// decomposition of the graph. The nodes in \c EVEN/D induce a subgraph kpeter@590: /// with factor-critical components, the nodes in \c ODD/A form the kpeter@590: /// canonical barrier, and the nodes in \c MATCHED/C induce a graph having kpeter@590: /// a perfect matching. The number of the factor-critical components deba@327: /// minus the number of barrier nodes is a lower bound on the alpar@330: /// unmatched nodes, and the matching is optimal if and only if this bound is kpeter@593: /// tight. This decomposition can be obtained using \ref status() or kpeter@593: /// \ref statusMap() after running the algorithm. deba@326: /// kpeter@593: /// \tparam GR The undirected graph type the algorithm runs on. kpeter@559: template deba@326: class MaxMatching { deba@327: public: deba@327: kpeter@590: /// The graph type of the algorithm kpeter@559: typedef GR Graph; kpeter@593: /// The type of the matching map deba@327: typedef typename Graph::template NodeMap deba@327: MatchingMap; deba@327: kpeter@590: ///\brief Status constants for Gallai-Edmonds decomposition. deba@327: /// deba@868: ///These constants are used for indicating the Gallai-Edmonds kpeter@590: ///decomposition of a graph. The nodes with status \c EVEN (or \c D) kpeter@590: ///induce a subgraph with factor-critical components, the nodes with kpeter@590: ///status \c ODD (or \c A) form the canonical barrier, and the nodes deba@868: ///with status \c MATCHED (or \c C) induce a subgraph having a kpeter@590: ///perfect matching. deba@327: enum Status { kpeter@590: EVEN = 1, ///< = 1. (\c D is an alias for \c EVEN.) kpeter@590: D = 1, kpeter@590: MATCHED = 0, ///< = 0. (\c C is an alias for \c MATCHED.) kpeter@590: C = 0, kpeter@590: ODD = -1, ///< = -1. (\c A is an alias for \c ODD.) kpeter@590: A = -1, kpeter@590: UNMATCHED = -2 ///< = -2. deba@327: }; deba@327: kpeter@593: /// The type of the status map deba@327: typedef typename Graph::template NodeMap StatusMap; deba@327: deba@327: private: deba@326: deba@326: TEMPLATE_GRAPH_TYPEDEFS(Graph); deba@326: deba@327: typedef UnionFindEnum BlossomSet; deba@327: typedef ExtendFindEnum TreeSet; deba@327: typedef RangeMap NodeIntMap; deba@327: typedef MatchingMap EarMap; deba@327: typedef std::vector NodeQueue; deba@327: deba@327: const Graph& _graph; deba@327: MatchingMap* _matching; deba@327: StatusMap* _status; deba@327: deba@327: EarMap* _ear; deba@327: deba@327: IntNodeMap* _blossom_set_index; deba@327: BlossomSet* _blossom_set; deba@327: NodeIntMap* _blossom_rep; deba@327: deba@327: IntNodeMap* _tree_set_index; deba@327: TreeSet* _tree_set; deba@327: deba@327: NodeQueue _node_queue; deba@327: int _process, _postpone, _last; deba@327: deba@1203: int _node_num, _unmatched; deba@327: deba@327: private: deba@327: deba@327: void createStructures() { deba@327: _node_num = countNodes(_graph); deba@327: if (!_matching) { deba@327: _matching = new MatchingMap(_graph); deba@327: } deba@327: if (!_status) { deba@327: _status = new StatusMap(_graph); deba@327: } deba@327: if (!_ear) { deba@327: _ear = new EarMap(_graph); deba@327: } deba@327: if (!_blossom_set) { deba@327: _blossom_set_index = new IntNodeMap(_graph); deba@327: _blossom_set = new BlossomSet(*_blossom_set_index); deba@327: } deba@327: if (!_blossom_rep) { deba@327: _blossom_rep = new NodeIntMap(_node_num); deba@327: } deba@327: if (!_tree_set) { deba@327: _tree_set_index = new IntNodeMap(_graph); deba@327: _tree_set = new TreeSet(*_tree_set_index); deba@327: } deba@327: _node_queue.resize(_node_num); deba@327: } deba@327: deba@327: void destroyStructures() { deba@327: if (_matching) { deba@327: delete _matching; deba@327: } deba@327: if (_status) { deba@327: delete _status; deba@327: } deba@327: if (_ear) { deba@327: delete _ear; deba@327: } deba@327: if (_blossom_set) { deba@327: delete _blossom_set; deba@327: delete _blossom_set_index; deba@327: } deba@327: if (_blossom_rep) { deba@327: delete _blossom_rep; deba@327: } deba@327: if (_tree_set) { deba@327: delete _tree_set_index; deba@327: delete _tree_set; deba@327: } deba@327: } deba@327: deba@327: void processDense(const Node& n) { deba@327: _process = _postpone = _last = 0; deba@327: _node_queue[_last++] = n; deba@327: deba@327: while (_process != _last) { deba@327: Node u = _node_queue[_process++]; deba@327: for (OutArcIt a(_graph, u); a != INVALID; ++a) { deba@327: Node v = _graph.target(a); deba@327: if ((*_status)[v] == MATCHED) { deba@327: extendOnArc(a); deba@327: } else if ((*_status)[v] == UNMATCHED) { deba@327: augmentOnArc(a); deba@1203: --_unmatched; deba@327: return; deba@327: } deba@327: } deba@327: } deba@327: deba@327: while (_postpone != _last) { deba@327: Node u = _node_queue[_postpone++]; deba@327: deba@327: for (OutArcIt a(_graph, u); a != INVALID ; ++a) { deba@327: Node v = _graph.target(a); deba@327: deba@327: if ((*_status)[v] == EVEN) { deba@327: if (_blossom_set->find(u) != _blossom_set->find(v)) { deba@327: shrinkOnEdge(a); deba@327: } deba@327: } deba@327: deba@327: while (_process != _last) { deba@327: Node w = _node_queue[_process++]; deba@327: for (OutArcIt b(_graph, w); b != INVALID; ++b) { deba@327: Node x = _graph.target(b); deba@327: if ((*_status)[x] == MATCHED) { deba@327: extendOnArc(b); deba@327: } else if ((*_status)[x] == UNMATCHED) { deba@327: augmentOnArc(b); deba@1203: --_unmatched; deba@327: return; deba@327: } deba@327: } deba@327: } deba@327: } deba@327: } deba@327: } deba@327: deba@327: void processSparse(const Node& n) { deba@327: _process = _last = 0; deba@327: _node_queue[_last++] = n; deba@327: while (_process != _last) { deba@327: Node u = _node_queue[_process++]; deba@327: for (OutArcIt a(_graph, u); a != INVALID; ++a) { deba@327: Node v = _graph.target(a); deba@327: deba@327: if ((*_status)[v] == EVEN) { deba@327: if (_blossom_set->find(u) != _blossom_set->find(v)) { deba@327: shrinkOnEdge(a); deba@327: } deba@327: } else if ((*_status)[v] == MATCHED) { deba@327: extendOnArc(a); deba@327: } else if ((*_status)[v] == UNMATCHED) { deba@327: augmentOnArc(a); deba@1203: --_unmatched; deba@327: return; deba@327: } deba@327: } deba@327: } deba@327: } deba@327: deba@327: void shrinkOnEdge(const Edge& e) { deba@327: Node nca = INVALID; deba@327: deba@327: { deba@327: std::set left_set, right_set; deba@327: deba@327: Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))]; deba@327: left_set.insert(left); deba@327: deba@327: Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))]; deba@327: right_set.insert(right); deba@327: deba@327: while (true) { deba@327: if ((*_matching)[left] == INVALID) break; deba@327: left = _graph.target((*_matching)[left]); deba@327: left = (*_blossom_rep)[_blossom_set-> deba@327: find(_graph.target((*_ear)[left]))]; deba@327: if (right_set.find(left) != right_set.end()) { deba@327: nca = left; deba@327: break; deba@327: } deba@327: left_set.insert(left); deba@327: deba@327: if ((*_matching)[right] == INVALID) break; deba@327: right = _graph.target((*_matching)[right]); deba@327: right = (*_blossom_rep)[_blossom_set-> deba@327: find(_graph.target((*_ear)[right]))]; deba@327: if (left_set.find(right) != left_set.end()) { deba@327: nca = right; deba@327: break; deba@327: } deba@327: right_set.insert(right); deba@327: } deba@327: deba@327: if (nca == INVALID) { deba@327: if ((*_matching)[left] == INVALID) { deba@327: nca = right; deba@327: while (left_set.find(nca) == left_set.end()) { deba@327: nca = _graph.target((*_matching)[nca]); deba@327: nca =(*_blossom_rep)[_blossom_set-> deba@327: find(_graph.target((*_ear)[nca]))]; deba@327: } deba@327: } else { deba@327: nca = left; deba@327: while (right_set.find(nca) == right_set.end()) { deba@327: nca = _graph.target((*_matching)[nca]); deba@327: nca = (*_blossom_rep)[_blossom_set-> deba@327: find(_graph.target((*_ear)[nca]))]; deba@327: } deba@327: } deba@327: } deba@327: } deba@327: deba@327: { deba@327: deba@327: Node node = _graph.u(e); deba@327: Arc arc = _graph.direct(e, true); deba@327: Node base = (*_blossom_rep)[_blossom_set->find(node)]; deba@327: deba@327: while (base != nca) { kpeter@581: (*_ear)[node] = arc; deba@327: deba@327: Node n = node; deba@327: while (n != base) { deba@327: n = _graph.target((*_matching)[n]); deba@327: Arc a = (*_ear)[n]; deba@327: n = _graph.target(a); kpeter@581: (*_ear)[n] = _graph.oppositeArc(a); deba@327: } deba@327: node = _graph.target((*_matching)[base]); deba@327: _tree_set->erase(base); deba@327: _tree_set->erase(node); deba@327: _blossom_set->insert(node, _blossom_set->find(base)); kpeter@581: (*_status)[node] = EVEN; deba@327: _node_queue[_last++] = node; deba@327: arc = _graph.oppositeArc((*_ear)[node]); deba@327: node = _graph.target((*_ear)[node]); deba@327: base = (*_blossom_rep)[_blossom_set->find(node)]; deba@327: _blossom_set->join(_graph.target(arc), base); deba@327: } deba@327: } deba@327: kpeter@581: (*_blossom_rep)[_blossom_set->find(nca)] = nca; deba@327: deba@327: { deba@327: deba@327: Node node = _graph.v(e); deba@327: Arc arc = _graph.direct(e, false); deba@327: Node base = (*_blossom_rep)[_blossom_set->find(node)]; deba@327: deba@327: while (base != nca) { kpeter@581: (*_ear)[node] = arc; deba@327: deba@327: Node n = node; deba@327: while (n != base) { deba@327: n = _graph.target((*_matching)[n]); deba@327: Arc a = (*_ear)[n]; deba@327: n = _graph.target(a); kpeter@581: (*_ear)[n] = _graph.oppositeArc(a); deba@327: } deba@327: node = _graph.target((*_matching)[base]); deba@327: _tree_set->erase(base); deba@327: _tree_set->erase(node); deba@327: _blossom_set->insert(node, _blossom_set->find(base)); kpeter@581: (*_status)[node] = EVEN; deba@327: _node_queue[_last++] = node; deba@327: arc = _graph.oppositeArc((*_ear)[node]); deba@327: node = _graph.target((*_ear)[node]); deba@327: base = (*_blossom_rep)[_blossom_set->find(node)]; deba@327: _blossom_set->join(_graph.target(arc), base); deba@327: } deba@327: } deba@327: kpeter@581: (*_blossom_rep)[_blossom_set->find(nca)] = nca; deba@327: } deba@327: deba@327: void extendOnArc(const Arc& a) { deba@327: Node base = _graph.source(a); deba@327: Node odd = _graph.target(a); deba@327: kpeter@581: (*_ear)[odd] = _graph.oppositeArc(a); deba@327: Node even = _graph.target((*_matching)[odd]); kpeter@581: (*_blossom_rep)[_blossom_set->insert(even)] = even; kpeter@581: (*_status)[odd] = ODD; kpeter@581: (*_status)[even] = EVEN; deba@327: int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]); deba@327: _tree_set->insert(odd, tree); deba@327: _tree_set->insert(even, tree); deba@327: _node_queue[_last++] = even; deba@327: deba@327: } deba@327: deba@327: void augmentOnArc(const Arc& a) { deba@327: Node even = _graph.source(a); deba@327: Node odd = _graph.target(a); deba@327: deba@327: int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]); deba@327: kpeter@581: (*_matching)[odd] = _graph.oppositeArc(a); kpeter@581: (*_status)[odd] = MATCHED; deba@327: deba@327: Arc arc = (*_matching)[even]; kpeter@581: (*_matching)[even] = a; deba@327: deba@327: while (arc != INVALID) { deba@327: odd = _graph.target(arc); deba@327: arc = (*_ear)[odd]; deba@327: even = _graph.target(arc); kpeter@581: (*_matching)[odd] = arc; deba@327: arc = (*_matching)[even]; kpeter@581: (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]); deba@327: } deba@327: deba@327: for (typename TreeSet::ItemIt it(*_tree_set, tree); deba@327: it != INVALID; ++it) { deba@327: if ((*_status)[it] == ODD) { kpeter@581: (*_status)[it] = MATCHED; deba@327: } else { deba@327: int blossom = _blossom_set->find(it); deba@327: for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom); deba@327: jt != INVALID; ++jt) { kpeter@581: (*_status)[jt] = MATCHED; deba@327: } deba@327: _blossom_set->eraseClass(blossom); deba@327: } deba@327: } deba@327: _tree_set->eraseClass(tree); deba@327: deba@327: } deba@326: deba@326: public: deba@326: deba@327: /// \brief Constructor deba@326: /// deba@327: /// Constructor. deba@327: MaxMatching(const Graph& graph) deba@327: : _graph(graph), _matching(0), _status(0), _ear(0), deba@327: _blossom_set_index(0), _blossom_set(0), _blossom_rep(0), deba@327: _tree_set_index(0), _tree_set(0) {} deba@327: deba@327: ~MaxMatching() { deba@327: destroyStructures(); deba@327: } deba@327: kpeter@590: /// \name Execution Control alpar@330: /// The simplest way to execute the algorithm is to use the kpeter@590: /// \c run() member function.\n kpeter@590: /// If you need better control on the execution, you have to call kpeter@590: /// one of the functions \ref init(), \ref greedyInit() or kpeter@590: /// \ref matchingInit() first, then you can start the algorithm with kpeter@590: /// \ref startSparse() or \ref startDense(). deba@327: deba@327: ///@{ deba@327: kpeter@590: /// \brief Set the initial matching to the empty matching. deba@326: /// kpeter@590: /// This function sets the initial matching to the empty matching. deba@326: void init() { deba@327: createStructures(); deba@327: for(NodeIt n(_graph); n != INVALID; ++n) { kpeter@581: (*_matching)[n] = INVALID; kpeter@581: (*_status)[n] = UNMATCHED; deba@326: } deba@1203: _unmatched = _node_num; deba@326: } deba@326: kpeter@590: /// \brief Find an initial matching in a greedy way. deba@326: /// kpeter@590: /// This function finds an initial matching in a greedy way. deba@326: void greedyInit() { deba@327: createStructures(); deba@327: for (NodeIt n(_graph); n != INVALID; ++n) { kpeter@581: (*_matching)[n] = INVALID; kpeter@581: (*_status)[n] = UNMATCHED; deba@326: } deba@1203: _unmatched = _node_num; deba@327: for (NodeIt n(_graph); n != INVALID; ++n) { deba@327: if ((*_matching)[n] == INVALID) { deba@327: for (OutArcIt a(_graph, n); a != INVALID ; ++a) { deba@327: Node v = _graph.target(a); deba@327: if ((*_matching)[v] == INVALID && v != n) { kpeter@581: (*_matching)[n] = a; kpeter@581: (*_status)[n] = MATCHED; kpeter@581: (*_matching)[v] = _graph.oppositeArc(a); kpeter@581: (*_status)[v] = MATCHED; deba@1203: _unmatched -= 2; deba@326: break; deba@326: } deba@326: } deba@326: } deba@326: } deba@326: } deba@326: deba@327: kpeter@590: /// \brief Initialize the matching from a map. deba@326: /// kpeter@590: /// This function initializes the matching from a \c bool valued edge kpeter@590: /// map. This map should have the property that there are no two incident kpeter@590: /// edges with \c true value, i.e. it really contains a matching. kpeter@559: /// \return \c true if the map contains a matching. deba@327: template deba@327: bool matchingInit(const MatchingMap& matching) { deba@327: createStructures(); deba@327: deba@327: for (NodeIt n(_graph); n != INVALID; ++n) { kpeter@581: (*_matching)[n] = INVALID; kpeter@581: (*_status)[n] = UNMATCHED; deba@326: } deba@1203: _unmatched = _node_num; deba@327: for(EdgeIt e(_graph); e!=INVALID; ++e) { deba@327: if (matching[e]) { deba@327: deba@327: Node u = _graph.u(e); deba@327: if ((*_matching)[u] != INVALID) return false; kpeter@581: (*_matching)[u] = _graph.direct(e, true); kpeter@581: (*_status)[u] = MATCHED; deba@327: deba@327: Node v = _graph.v(e); deba@327: if ((*_matching)[v] != INVALID) return false; kpeter@581: (*_matching)[v] = _graph.direct(e, false); kpeter@581: (*_status)[v] = MATCHED; deba@1203: deba@1203: _unmatched -= 2; deba@327: } deba@327: } deba@327: return true; deba@326: } deba@326: kpeter@590: /// \brief Start Edmonds' algorithm deba@326: /// deba@1203: /// This function runs the original Edmonds' algorithm. If the deba@1203: /// \c decomposition parameter is set to false, then the Gallai-Edmonds deba@1203: /// decomposition is not computed. kpeter@590: /// kpeter@651: /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be kpeter@590: /// called before using this function. deba@1203: void startSparse(bool decomposition = true) { deba@1203: int _unmatched_limit = decomposition ? 0 : 1; deba@1203: for(NodeIt n(_graph); _unmatched > _unmatched_limit; ++n) { deba@327: if ((*_status)[n] == UNMATCHED) { deba@327: (*_blossom_rep)[_blossom_set->insert(n)] = n; deba@327: _tree_set->insert(n); kpeter@581: (*_status)[n] = EVEN; deba@1203: --_unmatched; deba@327: processSparse(n); deba@326: } deba@326: } deba@326: } deba@326: deba@868: /// \brief Start Edmonds' algorithm with a heuristic improvement kpeter@590: /// for dense graphs deba@326: /// kpeter@590: /// This function runs Edmonds' algorithm with a heuristic of postponing deba@1203: /// shrinks, therefore resulting in a faster algorithm for dense graphs. If deba@1203: /// the \c decomposition parameter is set to false, then the Gallai-Edmonds deba@1203: /// decomposition is not computed. kpeter@590: /// kpeter@651: /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be kpeter@590: /// called before using this function. deba@1203: void startDense(bool decomposition = true) { deba@1203: int _unmatched_limit = decomposition ? 0 : 1; deba@1203: for(NodeIt n(_graph); _unmatched > _unmatched_limit; ++n) { deba@327: if ((*_status)[n] == UNMATCHED) { deba@327: (*_blossom_rep)[_blossom_set->insert(n)] = n; deba@327: _tree_set->insert(n); kpeter@581: (*_status)[n] = EVEN; deba@1203: --_unmatched; deba@327: processDense(n); deba@327: } deba@327: } deba@327: } deba@327: deba@327: kpeter@590: /// \brief Run Edmonds' algorithm deba@327: /// deba@868: /// This function runs Edmonds' algorithm. An additional heuristic of deba@1203: /// postponing shrinks is used for relatively dense graphs (for which deba@1203: /// m>=2*n holds). If the \c decomposition parameter is set to deba@1203: /// false, then the Gallai-Edmonds decomposition is not computed. In some deba@1203: /// cases, this can speed up the algorithm significantly, especially when a deba@1203: /// maximum matching is computed in a dense graph with odd number of nodes. deba@1203: void run(bool decomposition = true) { deba@327: if (countEdges(_graph) < 2 * countNodes(_graph)) { deba@326: greedyInit(); deba@1203: startSparse(decomposition); deba@326: } else { deba@326: init(); deba@1203: startDense(decomposition); deba@326: } deba@326: } deba@326: deba@327: /// @} deba@327: kpeter@590: /// \name Primal Solution kpeter@590: /// Functions to get the primal solution, i.e. the maximum matching. deba@327: deba@327: /// @{ deba@326: kpeter@590: /// \brief Return the size (cardinality) of the matching. deba@326: /// deba@868: /// This function returns the size (cardinality) of the current matching. kpeter@590: /// After run() it returns the size of the maximum matching in the graph. deba@327: int matchingSize() const { deba@327: int size = 0; deba@327: for (NodeIt n(_graph); n != INVALID; ++n) { deba@327: if ((*_matching)[n] != INVALID) { deba@327: ++size; deba@326: } deba@326: } deba@327: return size / 2; deba@326: } deba@326: kpeter@590: /// \brief Return \c true if the given edge is in the matching. deba@327: /// deba@868: /// This function returns \c true if the given edge is in the current kpeter@590: /// matching. deba@327: bool matching(const Edge& edge) const { deba@327: return edge == (*_matching)[_graph.u(edge)]; deba@327: } deba@327: kpeter@590: /// \brief Return the matching arc (or edge) incident to the given node. deba@327: /// kpeter@590: /// This function returns the matching arc (or edge) incident to the deba@868: /// given node in the current matching or \c INVALID if the node is kpeter@590: /// not covered by the matching. deba@327: Arc matching(const Node& n) const { deba@327: return (*_matching)[n]; deba@327: } deba@326: kpeter@593: /// \brief Return a const reference to the matching map. kpeter@593: /// kpeter@593: /// This function returns a const reference to a node map that stores kpeter@593: /// the matching arc (or edge) incident to each node. kpeter@593: const MatchingMap& matchingMap() const { kpeter@593: return *_matching; kpeter@593: } kpeter@593: kpeter@590: /// \brief Return the mate of the given node. deba@326: /// deba@868: /// This function returns the mate of the given node in the current kpeter@590: /// matching or \c INVALID if the node is not covered by the matching. deba@327: Node mate(const Node& n) const { deba@327: return (*_matching)[n] != INVALID ? deba@327: _graph.target((*_matching)[n]) : INVALID; deba@326: } deba@326: deba@327: /// @} deba@327: kpeter@590: /// \name Dual Solution deba@868: /// Functions to get the dual solution, i.e. the Gallai-Edmonds kpeter@590: /// decomposition. deba@327: deba@327: /// @{ deba@326: kpeter@590: /// \brief Return the status of the given node in the Edmonds-Gallai deba@326: /// decomposition. deba@326: /// kpeter@590: /// This function returns the \ref Status "status" of the given node kpeter@590: /// in the Edmonds-Gallai decomposition. kpeter@593: Status status(const Node& n) const { deba@327: return (*_status)[n]; deba@326: } deba@326: kpeter@593: /// \brief Return a const reference to the status map, which stores kpeter@593: /// the Edmonds-Gallai decomposition. kpeter@593: /// kpeter@593: /// This function returns a const reference to a node map that stores the kpeter@593: /// \ref Status "status" of each node in the Edmonds-Gallai decomposition. kpeter@593: const StatusMap& statusMap() const { kpeter@593: return *_status; kpeter@593: } kpeter@593: kpeter@590: /// \brief Return \c true if the given node is in the barrier. deba@326: /// kpeter@590: /// This function returns \c true if the given node is in the barrier. deba@327: bool barrier(const Node& n) const { deba@327: return (*_status)[n] == ODD; deba@326: } deba@326: deba@327: /// @} deba@326: deba@326: }; deba@326: deba@326: /// \ingroup matching deba@326: /// deba@326: /// \brief Weighted matching in general graphs deba@326: /// deba@326: /// This class provides an efficient implementation of Edmond's deba@326: /// maximum weighted matching algorithm. The implementation is based deba@326: /// on extensive use of priority queues and provides kpeter@559: /// \f$O(nm\log n)\f$ time complexity. deba@326: /// deba@868: /// The maximum weighted matching problem is to find a subset of the deba@868: /// edges in an undirected graph with maximum overall weight for which kpeter@590: /// each node has at most one incident edge. kpeter@590: /// It can be formulated with the following linear program. deba@326: /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f] deba@327: /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} deba@327: \quad \forall B\in\mathcal{O}\f] */ deba@326: /// \f[x_e \ge 0\quad \forall e\in E\f] deba@326: /// \f[\max \sum_{e\in E}x_ew_e\f] deba@327: /// where \f$\delta(X)\f$ is the set of edges incident to a node in deba@327: /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in deba@327: /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality deba@327: /// subsets of the nodes. deba@326: /// deba@326: /// The algorithm calculates an optimal matching and a proof of the deba@326: /// optimality. The solution of the dual problem can be used to check deba@327: /// the result of the algorithm. The dual linear problem is the kpeter@590: /// following. deba@327: /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)} deba@327: z_B \ge w_{uv} \quad \forall uv\in E\f] */ deba@326: /// \f[y_u \ge 0 \quad \forall u \in V\f] deba@326: /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] deba@327: /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} deba@327: \frac{\vert B \vert - 1}{2}z_B\f] */ deba@326: /// deba@868: /// The algorithm can be executed with the run() function. kpeter@590: /// After it the matching (the primal solution) and the dual solution deba@868: /// can be obtained using the query functions and the deba@868: /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, deba@868: /// which is able to iterate on the nodes of a blossom. kpeter@590: /// If the value type is integer, then the dual solution is multiplied kpeter@590: /// by \ref MaxWeightedMatching::dualScale "4". kpeter@590: /// kpeter@593: /// \tparam GR The undirected graph type the algorithm runs on. deba@868: /// \tparam WM The type edge weight map. The default type is kpeter@590: /// \ref concepts::Graph::EdgeMap "GR::EdgeMap". kpeter@590: #ifdef DOXYGEN kpeter@590: template kpeter@590: #else kpeter@559: template > kpeter@590: #endif deba@326: class MaxWeightedMatching { deba@326: public: deba@326: kpeter@590: /// The graph type of the algorithm kpeter@559: typedef GR Graph; kpeter@590: /// The type of the edge weight map kpeter@559: typedef WM WeightMap; kpeter@590: /// The value type of the edge weights deba@326: typedef typename WeightMap::Value Value; deba@326: kpeter@593: /// The type of the matching map kpeter@590: typedef typename Graph::template NodeMap kpeter@590: MatchingMap; kpeter@590: deba@326: /// \brief Scaling factor for dual solution deba@326: /// kpeter@590: /// Scaling factor for dual solution. It is equal to 4 or 1 deba@326: /// according to the value type. deba@326: static const int dualScale = deba@326: std::numeric_limits::is_integer ? 4 : 1; deba@326: deba@326: private: deba@326: deba@326: TEMPLATE_GRAPH_TYPEDEFS(Graph); deba@326: deba@326: typedef typename Graph::template NodeMap NodePotential; deba@326: typedef std::vector BlossomNodeList; deba@326: deba@326: struct BlossomVariable { deba@326: int begin, end; deba@326: Value value; deba@326: deba@326: BlossomVariable(int _begin, int _end, Value _value) deba@326: : begin(_begin), end(_end), value(_value) {} deba@326: deba@326: }; deba@326: deba@326: typedef std::vector BlossomPotential; deba@326: deba@326: const Graph& _graph; deba@326: const WeightMap& _weight; deba@326: deba@326: MatchingMap* _matching; deba@326: deba@326: NodePotential* _node_potential; deba@326: deba@326: BlossomPotential _blossom_potential; deba@326: BlossomNodeList _blossom_node_list; deba@326: deba@326: int _node_num; deba@326: int _blossom_num; deba@326: deba@326: typedef RangeMap IntIntMap; deba@326: deba@326: enum Status { deba@868: EVEN = -1, MATCHED = 0, ODD = 1 deba@326: }; deba@326: deba@327: typedef HeapUnionFind BlossomSet; deba@326: struct BlossomData { deba@326: int tree; deba@326: Status status; deba@326: Arc pred, next; deba@326: Value pot, offset; deba@326: Node base; deba@326: }; deba@326: deba@327: IntNodeMap *_blossom_index; deba@326: BlossomSet *_blossom_set; deba@326: RangeMap* _blossom_data; deba@326: deba@327: IntNodeMap *_node_index; deba@327: IntArcMap *_node_heap_index; deba@326: deba@326: struct NodeData { deba@326: deba@327: NodeData(IntArcMap& node_heap_index) deba@326: : heap(node_heap_index) {} deba@326: deba@326: int blossom; deba@326: Value pot; deba@327: BinHeap heap; deba@326: std::map heap_index; deba@326: deba@326: int tree; deba@326: }; deba@326: deba@326: RangeMap* _node_data; deba@326: deba@326: typedef ExtendFindEnum TreeSet; deba@326: deba@326: IntIntMap *_tree_set_index; deba@326: TreeSet *_tree_set; deba@326: deba@327: IntNodeMap *_delta1_index; deba@327: BinHeap *_delta1; deba@326: deba@326: IntIntMap *_delta2_index; deba@326: BinHeap *_delta2; deba@326: deba@327: IntEdgeMap *_delta3_index; deba@327: BinHeap *_delta3; deba@326: deba@326: IntIntMap *_delta4_index; deba@326: BinHeap *_delta4; deba@326: deba@326: Value _delta_sum; deba@870: int _unmatched; deba@870: deba@870: typedef MaxWeightedFractionalMatching FractionalMatching; deba@870: FractionalMatching *_fractional; deba@326: deba@326: void createStructures() { deba@326: _node_num = countNodes(_graph); deba@326: _blossom_num = _node_num * 3 / 2; deba@326: deba@326: if (!_matching) { deba@326: _matching = new MatchingMap(_graph); deba@326: } deba@867: deba@326: if (!_node_potential) { deba@326: _node_potential = new NodePotential(_graph); deba@326: } deba@867: deba@326: if (!_blossom_set) { deba@327: _blossom_index = new IntNodeMap(_graph); deba@326: _blossom_set = new BlossomSet(*_blossom_index); deba@326: _blossom_data = new RangeMap(_blossom_num); deba@867: } else if (_blossom_data->size() != _blossom_num) { deba@867: delete _blossom_data; deba@867: _blossom_data = new RangeMap(_blossom_num); deba@326: } deba@326: deba@326: if (!_node_index) { deba@327: _node_index = new IntNodeMap(_graph); deba@327: _node_heap_index = new IntArcMap(_graph); deba@326: _node_data = new RangeMap(_node_num, deba@867: NodeData(*_node_heap_index)); deba@867: } else { deba@867: delete _node_data; deba@867: _node_data = new RangeMap(_node_num, deba@867: NodeData(*_node_heap_index)); deba@326: } deba@326: deba@326: if (!_tree_set) { deba@326: _tree_set_index = new IntIntMap(_blossom_num); deba@326: _tree_set = new TreeSet(*_tree_set_index); deba@867: } else { deba@867: _tree_set_index->resize(_blossom_num); deba@326: } deba@867: deba@326: if (!_delta1) { deba@327: _delta1_index = new IntNodeMap(_graph); deba@327: _delta1 = new BinHeap(*_delta1_index); deba@326: } deba@867: deba@326: if (!_delta2) { deba@326: _delta2_index = new IntIntMap(_blossom_num); deba@326: _delta2 = new BinHeap(*_delta2_index); deba@867: } else { deba@867: _delta2_index->resize(_blossom_num); deba@326: } deba@867: deba@326: if (!_delta3) { deba@327: _delta3_index = new IntEdgeMap(_graph); deba@327: _delta3 = new BinHeap(*_delta3_index); deba@326: } deba@867: deba@326: if (!_delta4) { deba@326: _delta4_index = new IntIntMap(_blossom_num); deba@326: _delta4 = new BinHeap(*_delta4_index); deba@867: } else { deba@867: _delta4_index->resize(_blossom_num); deba@326: } deba@326: } deba@326: deba@326: void destroyStructures() { deba@326: if (_matching) { deba@326: delete _matching; deba@326: } deba@326: if (_node_potential) { deba@326: delete _node_potential; deba@326: } deba@326: if (_blossom_set) { deba@326: delete _blossom_index; deba@326: delete _blossom_set; deba@326: delete _blossom_data; deba@326: } deba@326: deba@326: if (_node_index) { deba@326: delete _node_index; deba@326: delete _node_heap_index; deba@326: delete _node_data; deba@326: } deba@326: deba@326: if (_tree_set) { deba@326: delete _tree_set_index; deba@326: delete _tree_set; deba@326: } deba@326: if (_delta1) { deba@326: delete _delta1_index; deba@326: delete _delta1; deba@326: } deba@326: if (_delta2) { deba@326: delete _delta2_index; deba@326: delete _delta2; deba@326: } deba@326: if (_delta3) { deba@326: delete _delta3_index; deba@326: delete _delta3; deba@326: } deba@326: if (_delta4) { deba@326: delete _delta4_index; deba@326: delete _delta4; deba@326: } deba@326: } deba@326: deba@326: void matchedToEven(int blossom, int tree) { deba@326: if (_delta2->state(blossom) == _delta2->IN_HEAP) { deba@326: _delta2->erase(blossom); deba@326: } deba@326: deba@326: if (!_blossom_set->trivial(blossom)) { deba@326: (*_blossom_data)[blossom].pot -= deba@326: 2 * (_delta_sum - (*_blossom_data)[blossom].offset); deba@326: } deba@326: deba@326: for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); deba@326: n != INVALID; ++n) { deba@326: deba@326: _blossom_set->increase(n, std::numeric_limits::max()); deba@326: int ni = (*_node_index)[n]; deba@326: deba@326: (*_node_data)[ni].heap.clear(); deba@326: (*_node_data)[ni].heap_index.clear(); deba@326: deba@326: (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset; deba@326: deba@326: _delta1->push(n, (*_node_data)[ni].pot); deba@326: deba@326: for (InArcIt e(_graph, n); e != INVALID; ++e) { deba@326: Node v = _graph.source(e); deba@326: int vb = _blossom_set->find(v); deba@326: int vi = (*_node_index)[v]; deba@326: deba@326: Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - deba@326: dualScale * _weight[e]; deba@326: deba@326: if ((*_blossom_data)[vb].status == EVEN) { deba@326: if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { deba@326: _delta3->push(e, rw / 2); deba@326: } deba@326: } else { deba@326: typename std::map::iterator it = deba@326: (*_node_data)[vi].heap_index.find(tree); deba@326: deba@326: if (it != (*_node_data)[vi].heap_index.end()) { deba@326: if ((*_node_data)[vi].heap[it->second] > rw) { deba@326: (*_node_data)[vi].heap.replace(it->second, e); deba@326: (*_node_data)[vi].heap.decrease(e, rw); deba@326: it->second = e; deba@326: } deba@326: } else { deba@326: (*_node_data)[vi].heap.push(e, rw); deba@326: (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); deba@326: } deba@326: deba@326: if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { deba@326: _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); deba@326: deba@326: if ((*_blossom_data)[vb].status == MATCHED) { deba@326: if (_delta2->state(vb) != _delta2->IN_HEAP) { deba@326: _delta2->push(vb, _blossom_set->classPrio(vb) - deba@326: (*_blossom_data)[vb].offset); deba@326: } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - deba@326: (*_blossom_data)[vb].offset) { deba@326: _delta2->decrease(vb, _blossom_set->classPrio(vb) - deba@326: (*_blossom_data)[vb].offset); deba@326: } deba@326: } deba@326: } deba@326: } deba@326: } deba@326: } deba@326: (*_blossom_data)[blossom].offset = 0; deba@326: } deba@326: deba@868: void matchedToOdd(int blossom) { deba@326: if (_delta2->state(blossom) == _delta2->IN_HEAP) { deba@326: _delta2->erase(blossom); deba@326: } deba@868: (*_blossom_data)[blossom].offset += _delta_sum; deba@868: if (!_blossom_set->trivial(blossom)) { deba@868: _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + deba@868: (*_blossom_data)[blossom].offset); deba@868: } deba@868: } deba@868: deba@868: void evenToMatched(int blossom, int tree) { deba@868: if (!_blossom_set->trivial(blossom)) { deba@868: (*_blossom_data)[blossom].pot += 2 * _delta_sum; deba@868: } deba@326: deba@326: for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); deba@326: n != INVALID; ++n) { deba@326: int ni = (*_node_index)[n]; deba@868: (*_node_data)[ni].pot -= _delta_sum; deba@868: deba@868: _delta1->erase(n); deba@868: deba@868: for (InArcIt e(_graph, n); e != INVALID; ++e) { deba@868: Node v = _graph.source(e); deba@326: int vb = _blossom_set->find(v); deba@326: int vi = (*_node_index)[v]; deba@326: deba@326: Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - deba@326: dualScale * _weight[e]; deba@326: deba@868: if (vb == blossom) { deba@868: if (_delta3->state(e) == _delta3->IN_HEAP) { deba@868: _delta3->erase(e); deba@868: } deba@868: } else if ((*_blossom_data)[vb].status == EVEN) { deba@868: deba@868: if (_delta3->state(e) == _delta3->IN_HEAP) { deba@868: _delta3->erase(e); deba@868: } deba@868: deba@868: int vt = _tree_set->find(vb); deba@868: deba@868: if (vt != tree) { deba@868: deba@868: Arc r = _graph.oppositeArc(e); deba@868: deba@868: typename std::map::iterator it = deba@868: (*_node_data)[ni].heap_index.find(vt); deba@868: deba@868: if (it != (*_node_data)[ni].heap_index.end()) { deba@868: if ((*_node_data)[ni].heap[it->second] > rw) { deba@868: (*_node_data)[ni].heap.replace(it->second, r); deba@868: (*_node_data)[ni].heap.decrease(r, rw); deba@868: it->second = r; deba@868: } deba@868: } else { deba@868: (*_node_data)[ni].heap.push(r, rw); deba@868: (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); deba@868: } deba@868: deba@868: if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { deba@868: _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); deba@868: deba@868: if (_delta2->state(blossom) != _delta2->IN_HEAP) { deba@868: _delta2->push(blossom, _blossom_set->classPrio(blossom) - deba@868: (*_blossom_data)[blossom].offset); deba@868: } else if ((*_delta2)[blossom] > deba@868: _blossom_set->classPrio(blossom) - deba@868: (*_blossom_data)[blossom].offset){ deba@868: _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - deba@868: (*_blossom_data)[blossom].offset); deba@868: } deba@868: } deba@868: } deba@868: } else { deba@868: deba@868: typename std::map::iterator it = deba@868: (*_node_data)[vi].heap_index.find(tree); deba@868: deba@868: if (it != (*_node_data)[vi].heap_index.end()) { deba@868: (*_node_data)[vi].heap.erase(it->second); deba@868: (*_node_data)[vi].heap_index.erase(it); deba@868: if ((*_node_data)[vi].heap.empty()) { deba@868: _blossom_set->increase(v, std::numeric_limits::max()); deba@868: } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { deba@868: _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); deba@868: } deba@868: deba@868: if ((*_blossom_data)[vb].status == MATCHED) { deba@868: if (_blossom_set->classPrio(vb) == deba@868: std::numeric_limits::max()) { deba@868: _delta2->erase(vb); deba@868: } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - deba@868: (*_blossom_data)[vb].offset) { deba@868: _delta2->increase(vb, _blossom_set->classPrio(vb) - deba@868: (*_blossom_data)[vb].offset); deba@868: } deba@868: } deba@326: } deba@326: } deba@326: } deba@326: } deba@326: } deba@326: deba@868: void oddToMatched(int blossom) { deba@868: (*_blossom_data)[blossom].offset -= _delta_sum; deba@868: deba@868: if (_blossom_set->classPrio(blossom) != deba@868: std::numeric_limits::max()) { deba@868: _delta2->push(blossom, _blossom_set->classPrio(blossom) - deba@868: (*_blossom_data)[blossom].offset); deba@868: } deba@868: deba@868: if (!_blossom_set->trivial(blossom)) { deba@868: _delta4->erase(blossom); deba@868: } deba@868: } deba@868: deba@868: void oddToEven(int blossom, int tree) { deba@868: if (!_blossom_set->trivial(blossom)) { deba@868: _delta4->erase(blossom); deba@868: (*_blossom_data)[blossom].pot -= deba@868: 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); deba@868: } deba@868: deba@326: for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); deba@326: n != INVALID; ++n) { deba@326: int ni = (*_node_index)[n]; deba@326: deba@868: _blossom_set->increase(n, std::numeric_limits::max()); deba@868: deba@868: (*_node_data)[ni].heap.clear(); deba@868: (*_node_data)[ni].heap_index.clear(); deba@868: (*_node_data)[ni].pot += deba@868: 2 * _delta_sum - (*_blossom_data)[blossom].offset; deba@868: deba@868: _delta1->push(n, (*_node_data)[ni].pot); deba@868: deba@326: for (InArcIt e(_graph, n); e != INVALID; ++e) { deba@326: Node v = _graph.source(e); deba@326: int vb = _blossom_set->find(v); deba@326: int vi = (*_node_index)[v]; deba@326: deba@326: Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - deba@326: dualScale * _weight[e]; deba@326: deba@868: if ((*_blossom_data)[vb].status == EVEN) { deba@868: if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { deba@868: _delta3->push(e, rw / 2); deba@326: } deba@868: } else { deba@326: deba@326: typename std::map::iterator it = deba@868: (*_node_data)[vi].heap_index.find(tree); deba@868: deba@868: if (it != (*_node_data)[vi].heap_index.end()) { deba@868: if ((*_node_data)[vi].heap[it->second] > rw) { deba@868: (*_node_data)[vi].heap.replace(it->second, e); deba@868: (*_node_data)[vi].heap.decrease(e, rw); deba@868: it->second = e; deba@326: } deba@326: } else { deba@868: (*_node_data)[vi].heap.push(e, rw); deba@868: (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); deba@326: } deba@326: deba@868: if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { deba@868: _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); deba@868: deba@868: if ((*_blossom_data)[vb].status == MATCHED) { deba@868: if (_delta2->state(vb) != _delta2->IN_HEAP) { deba@868: _delta2->push(vb, _blossom_set->classPrio(vb) - deba@868: (*_blossom_data)[vb].offset); deba@868: } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - deba@868: (*_blossom_data)[vb].offset) { deba@868: _delta2->decrease(vb, _blossom_set->classPrio(vb) - deba@868: (*_blossom_data)[vb].offset); deba@868: } deba@326: } deba@326: } deba@326: } deba@326: } deba@326: } deba@868: (*_blossom_data)[blossom].offset = 0; deba@326: } deba@326: deba@326: void alternatePath(int even, int tree) { deba@326: int odd; deba@326: deba@326: evenToMatched(even, tree); deba@326: (*_blossom_data)[even].status = MATCHED; deba@326: deba@326: while ((*_blossom_data)[even].pred != INVALID) { deba@326: odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred)); deba@326: (*_blossom_data)[odd].status = MATCHED; deba@326: oddToMatched(odd); deba@326: (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; deba@326: deba@326: even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred)); deba@326: (*_blossom_data)[even].status = MATCHED; deba@326: evenToMatched(even, tree); deba@326: (*_blossom_data)[even].next = deba@326: _graph.oppositeArc((*_blossom_data)[odd].pred); deba@326: } deba@326: deba@326: } deba@326: deba@326: void destroyTree(int tree) { deba@326: for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { deba@326: if ((*_blossom_data)[b].status == EVEN) { deba@326: (*_blossom_data)[b].status = MATCHED; deba@326: evenToMatched(b, tree); deba@326: } else if ((*_blossom_data)[b].status == ODD) { deba@326: (*_blossom_data)[b].status = MATCHED; deba@326: oddToMatched(b); deba@326: } deba@326: } deba@326: _tree_set->eraseClass(tree); deba@326: } deba@326: deba@326: deba@326: void unmatchNode(const Node& node) { deba@326: int blossom = _blossom_set->find(node); deba@326: int tree = _tree_set->find(blossom); deba@326: deba@326: alternatePath(blossom, tree); deba@326: destroyTree(tree); deba@326: deba@326: (*_blossom_data)[blossom].base = node; deba@868: (*_blossom_data)[blossom].next = INVALID; deba@326: } deba@326: deba@327: void augmentOnEdge(const Edge& edge) { deba@327: deba@327: int left = _blossom_set->find(_graph.u(edge)); deba@327: int right = _blossom_set->find(_graph.v(edge)); deba@326: deba@868: int left_tree = _tree_set->find(left); deba@868: alternatePath(left, left_tree); deba@868: destroyTree(left_tree); deba@868: deba@868: int right_tree = _tree_set->find(right); deba@868: alternatePath(right, right_tree); deba@868: destroyTree(right_tree); deba@326: deba@327: (*_blossom_data)[left].next = _graph.direct(edge, true); deba@327: (*_blossom_data)[right].next = _graph.direct(edge, false); deba@326: } deba@326: deba@868: void augmentOnArc(const Arc& arc) { deba@868: deba@868: int left = _blossom_set->find(_graph.source(arc)); deba@868: int right = _blossom_set->find(_graph.target(arc)); deba@868: deba@868: (*_blossom_data)[left].status = MATCHED; deba@868: deba@868: int right_tree = _tree_set->find(right); deba@868: alternatePath(right, right_tree); deba@868: destroyTree(right_tree); deba@868: deba@868: (*_blossom_data)[left].next = arc; deba@868: (*_blossom_data)[right].next = _graph.oppositeArc(arc); deba@868: } deba@868: deba@326: void extendOnArc(const Arc& arc) { deba@326: int base = _blossom_set->find(_graph.target(arc)); deba@326: int tree = _tree_set->find(base); deba@326: deba@326: int odd = _blossom_set->find(_graph.source(arc)); deba@326: _tree_set->insert(odd, tree); deba@326: (*_blossom_data)[odd].status = ODD; deba@326: matchedToOdd(odd); deba@326: (*_blossom_data)[odd].pred = arc; deba@326: deba@326: int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next)); deba@326: (*_blossom_data)[even].pred = (*_blossom_data)[even].next; deba@326: _tree_set->insert(even, tree); deba@326: (*_blossom_data)[even].status = EVEN; deba@326: matchedToEven(even, tree); deba@326: } deba@326: deba@327: void shrinkOnEdge(const Edge& edge, int tree) { deba@326: int nca = -1; deba@326: std::vector left_path, right_path; deba@326: deba@326: { deba@326: std::set left_set, right_set; deba@326: int left = _blossom_set->find(_graph.u(edge)); deba@326: left_path.push_back(left); deba@326: left_set.insert(left); deba@326: deba@326: int right = _blossom_set->find(_graph.v(edge)); deba@326: right_path.push_back(right); deba@326: right_set.insert(right); deba@326: deba@326: while (true) { deba@326: deba@326: if ((*_blossom_data)[left].pred == INVALID) break; deba@326: deba@326: left = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); deba@326: left_path.push_back(left); deba@326: left = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); deba@326: left_path.push_back(left); deba@326: deba@326: left_set.insert(left); deba@326: deba@326: if (right_set.find(left) != right_set.end()) { deba@326: nca = left; deba@326: break; deba@326: } deba@326: deba@326: if ((*_blossom_data)[right].pred == INVALID) break; deba@326: deba@326: right = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); deba@326: right_path.push_back(right); deba@326: right = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); deba@326: right_path.push_back(right); deba@326: deba@326: right_set.insert(right); deba@326: deba@326: if (left_set.find(right) != left_set.end()) { deba@326: nca = right; deba@326: break; deba@326: } deba@326: deba@326: } deba@326: deba@326: if (nca == -1) { deba@326: if ((*_blossom_data)[left].pred == INVALID) { deba@326: nca = right; deba@326: while (left_set.find(nca) == left_set.end()) { deba@326: nca = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); deba@326: right_path.push_back(nca); deba@326: nca = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); deba@326: right_path.push_back(nca); deba@326: } deba@326: } else { deba@326: nca = left; deba@326: while (right_set.find(nca) == right_set.end()) { deba@326: nca = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); deba@326: left_path.push_back(nca); deba@326: nca = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); deba@326: left_path.push_back(nca); deba@326: } deba@326: } deba@326: } deba@326: } deba@326: deba@326: std::vector subblossoms; deba@326: Arc prev; deba@326: deba@326: prev = _graph.direct(edge, true); deba@326: for (int i = 0; left_path[i] != nca; i += 2) { deba@326: subblossoms.push_back(left_path[i]); deba@326: (*_blossom_data)[left_path[i]].next = prev; deba@326: _tree_set->erase(left_path[i]); deba@326: deba@326: subblossoms.push_back(left_path[i + 1]); deba@326: (*_blossom_data)[left_path[i + 1]].status = EVEN; deba@326: oddToEven(left_path[i + 1], tree); deba@326: _tree_set->erase(left_path[i + 1]); deba@326: prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred); deba@326: } deba@326: deba@326: int k = 0; deba@326: while (right_path[k] != nca) ++k; deba@326: deba@326: subblossoms.push_back(nca); deba@326: (*_blossom_data)[nca].next = prev; deba@326: deba@326: for (int i = k - 2; i >= 0; i -= 2) { deba@326: subblossoms.push_back(right_path[i + 1]); deba@326: (*_blossom_data)[right_path[i + 1]].status = EVEN; deba@326: oddToEven(right_path[i + 1], tree); deba@326: _tree_set->erase(right_path[i + 1]); deba@326: deba@326: (*_blossom_data)[right_path[i + 1]].next = deba@326: (*_blossom_data)[right_path[i + 1]].pred; deba@326: deba@326: subblossoms.push_back(right_path[i]); deba@326: _tree_set->erase(right_path[i]); deba@326: } deba@326: deba@326: int surface = deba@326: _blossom_set->join(subblossoms.begin(), subblossoms.end()); deba@326: deba@326: for (int i = 0; i < int(subblossoms.size()); ++i) { deba@326: if (!_blossom_set->trivial(subblossoms[i])) { deba@326: (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; deba@326: } deba@326: (*_blossom_data)[subblossoms[i]].status = MATCHED; deba@326: } deba@326: deba@326: (*_blossom_data)[surface].pot = -2 * _delta_sum; deba@326: (*_blossom_data)[surface].offset = 0; deba@326: (*_blossom_data)[surface].status = EVEN; deba@326: (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; deba@326: (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; deba@326: deba@326: _tree_set->insert(surface, tree); deba@326: _tree_set->erase(nca); deba@326: } deba@326: deba@326: void splitBlossom(int blossom) { deba@326: Arc next = (*_blossom_data)[blossom].next; deba@326: Arc pred = (*_blossom_data)[blossom].pred; deba@326: deba@326: int tree = _tree_set->find(blossom); deba@326: deba@326: (*_blossom_data)[blossom].status = MATCHED; deba@326: oddToMatched(blossom); deba@326: if (_delta2->state(blossom) == _delta2->IN_HEAP) { deba@326: _delta2->erase(blossom); deba@326: } deba@326: deba@326: std::vector subblossoms; deba@326: _blossom_set->split(blossom, std::back_inserter(subblossoms)); deba@326: deba@326: Value offset = (*_blossom_data)[blossom].offset; deba@326: int b = _blossom_set->find(_graph.source(pred)); deba@326: int d = _blossom_set->find(_graph.source(next)); deba@326: deba@326: int ib = -1, id = -1; deba@326: for (int i = 0; i < int(subblossoms.size()); ++i) { deba@326: if (subblossoms[i] == b) ib = i; deba@326: if (subblossoms[i] == d) id = i; deba@326: deba@326: (*_blossom_data)[subblossoms[i]].offset = offset; deba@326: if (!_blossom_set->trivial(subblossoms[i])) { deba@326: (*_blossom_data)[subblossoms[i]].pot -= 2 * offset; deba@326: } deba@326: if (_blossom_set->classPrio(subblossoms[i]) != deba@326: std::numeric_limits::max()) { deba@326: _delta2->push(subblossoms[i], deba@326: _blossom_set->classPrio(subblossoms[i]) - deba@326: (*_blossom_data)[subblossoms[i]].offset); deba@326: } deba@326: } deba@326: deba@326: if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) { deba@326: for (int i = (id + 1) % subblossoms.size(); deba@326: i != ib; i = (i + 2) % subblossoms.size()) { deba@326: int sb = subblossoms[i]; deba@326: int tb = subblossoms[(i + 1) % subblossoms.size()]; deba@326: (*_blossom_data)[sb].next = deba@326: _graph.oppositeArc((*_blossom_data)[tb].next); deba@326: } deba@326: deba@326: for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { deba@326: int sb = subblossoms[i]; deba@326: int tb = subblossoms[(i + 1) % subblossoms.size()]; deba@326: int ub = subblossoms[(i + 2) % subblossoms.size()]; deba@326: deba@326: (*_blossom_data)[sb].status = ODD; deba@326: matchedToOdd(sb); deba@326: _tree_set->insert(sb, tree); deba@326: (*_blossom_data)[sb].pred = pred; deba@326: (*_blossom_data)[sb].next = deba@868: _graph.oppositeArc((*_blossom_data)[tb].next); deba@326: deba@326: pred = (*_blossom_data)[ub].next; deba@326: deba@326: (*_blossom_data)[tb].status = EVEN; deba@326: matchedToEven(tb, tree); deba@326: _tree_set->insert(tb, tree); deba@326: (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; deba@326: } deba@326: deba@326: (*_blossom_data)[subblossoms[id]].status = ODD; deba@326: matchedToOdd(subblossoms[id]); deba@326: _tree_set->insert(subblossoms[id], tree); deba@326: (*_blossom_data)[subblossoms[id]].next = next; deba@326: (*_blossom_data)[subblossoms[id]].pred = pred; deba@326: deba@326: } else { deba@326: deba@326: for (int i = (ib + 1) % subblossoms.size(); deba@326: i != id; i = (i + 2) % subblossoms.size()) { deba@326: int sb = subblossoms[i]; deba@326: int tb = subblossoms[(i + 1) % subblossoms.size()]; deba@326: (*_blossom_data)[sb].next = deba@326: _graph.oppositeArc((*_blossom_data)[tb].next); deba@326: } deba@326: deba@326: for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { deba@326: int sb = subblossoms[i]; deba@326: int tb = subblossoms[(i + 1) % subblossoms.size()]; deba@326: int ub = subblossoms[(i + 2) % subblossoms.size()]; deba@326: deba@326: (*_blossom_data)[sb].status = ODD; deba@326: matchedToOdd(sb); deba@326: _tree_set->insert(sb, tree); deba@326: (*_blossom_data)[sb].next = next; deba@326: (*_blossom_data)[sb].pred = deba@326: _graph.oppositeArc((*_blossom_data)[tb].next); deba@326: deba@326: (*_blossom_data)[tb].status = EVEN; deba@326: matchedToEven(tb, tree); deba@326: _tree_set->insert(tb, tree); deba@326: (*_blossom_data)[tb].pred = deba@326: (*_blossom_data)[tb].next = deba@326: _graph.oppositeArc((*_blossom_data)[ub].next); deba@326: next = (*_blossom_data)[ub].next; deba@326: } deba@326: deba@326: (*_blossom_data)[subblossoms[ib]].status = ODD; deba@326: matchedToOdd(subblossoms[ib]); deba@326: _tree_set->insert(subblossoms[ib], tree); deba@326: (*_blossom_data)[subblossoms[ib]].next = next; deba@326: (*_blossom_data)[subblossoms[ib]].pred = pred; deba@326: } deba@326: _tree_set->erase(blossom); deba@326: } deba@326: deba@326: void extractBlossom(int blossom, const Node& base, const Arc& matching) { deba@326: if (_blossom_set->trivial(blossom)) { deba@326: int bi = (*_node_index)[base]; deba@326: Value pot = (*_node_data)[bi].pot; deba@326: kpeter@581: (*_matching)[base] = matching; deba@326: _blossom_node_list.push_back(base); kpeter@581: (*_node_potential)[base] = pot; deba@326: } else { deba@326: deba@326: Value pot = (*_blossom_data)[blossom].pot; deba@326: int bn = _blossom_node_list.size(); deba@326: deba@326: std::vector subblossoms; deba@326: _blossom_set->split(blossom, std::back_inserter(subblossoms)); deba@326: int b = _blossom_set->find(base); deba@326: int ib = -1; deba@326: for (int i = 0; i < int(subblossoms.size()); ++i) { deba@326: if (subblossoms[i] == b) { ib = i; break; } deba@326: } deba@326: deba@326: for (int i = 1; i < int(subblossoms.size()); i += 2) { deba@326: int sb = subblossoms[(ib + i) % subblossoms.size()]; deba@326: int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; deba@326: deba@326: Arc m = (*_blossom_data)[tb].next; deba@326: extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m)); deba@326: extractBlossom(tb, _graph.source(m), m); deba@326: } deba@326: extractBlossom(subblossoms[ib], base, matching); deba@326: deba@326: int en = _blossom_node_list.size(); deba@326: deba@326: _blossom_potential.push_back(BlossomVariable(bn, en, pot)); deba@326: } deba@326: } deba@326: deba@326: void extractMatching() { deba@326: std::vector blossoms; deba@326: for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { deba@326: blossoms.push_back(c); deba@326: } deba@326: deba@326: for (int i = 0; i < int(blossoms.size()); ++i) { deba@868: if ((*_blossom_data)[blossoms[i]].next != INVALID) { deba@326: deba@326: Value offset = (*_blossom_data)[blossoms[i]].offset; deba@326: (*_blossom_data)[blossoms[i]].pot += 2 * offset; deba@326: for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); deba@326: n != INVALID; ++n) { deba@326: (*_node_data)[(*_node_index)[n]].pot -= offset; deba@326: } deba@326: deba@326: Arc matching = (*_blossom_data)[blossoms[i]].next; deba@326: Node base = _graph.source(matching); deba@326: extractBlossom(blossoms[i], base, matching); deba@326: } else { deba@326: Node base = (*_blossom_data)[blossoms[i]].base; deba@326: extractBlossom(blossoms[i], base, INVALID); deba@326: } deba@326: } deba@326: } deba@326: deba@326: public: deba@326: deba@326: /// \brief Constructor deba@326: /// deba@326: /// Constructor. deba@326: MaxWeightedMatching(const Graph& graph, const WeightMap& weight) deba@326: : _graph(graph), _weight(weight), _matching(0), deba@326: _node_potential(0), _blossom_potential(), _blossom_node_list(), deba@326: _node_num(0), _blossom_num(0), deba@326: deba@326: _blossom_index(0), _blossom_set(0), _blossom_data(0), deba@326: _node_index(0), _node_heap_index(0), _node_data(0), deba@326: _tree_set_index(0), _tree_set(0), deba@326: deba@326: _delta1_index(0), _delta1(0), deba@326: _delta2_index(0), _delta2(0), deba@326: _delta3_index(0), _delta3(0), deba@326: _delta4_index(0), _delta4(0), deba@326: deba@870: _delta_sum(), _unmatched(0), deba@870: deba@870: _fractional(0) deba@870: {} deba@326: deba@326: ~MaxWeightedMatching() { deba@326: destroyStructures(); deba@870: if (_fractional) { deba@870: delete _fractional; deba@870: } deba@326: } deba@326: kpeter@590: /// \name Execution Control alpar@330: /// The simplest way to execute the algorithm is to use the kpeter@590: /// \ref run() member function. deba@326: deba@326: ///@{ deba@326: deba@326: /// \brief Initialize the algorithm deba@326: /// kpeter@590: /// This function initializes the algorithm. deba@326: void init() { deba@326: createStructures(); deba@326: deba@867: _blossom_node_list.clear(); deba@867: _blossom_potential.clear(); deba@867: deba@326: for (ArcIt e(_graph); e != INVALID; ++e) { kpeter@581: (*_node_heap_index)[e] = BinHeap::PRE_HEAP; deba@326: } deba@326: for (NodeIt n(_graph); n != INVALID; ++n) { kpeter@581: (*_delta1_index)[n] = _delta1->PRE_HEAP; deba@326: } deba@326: for (EdgeIt e(_graph); e != INVALID; ++e) { kpeter@581: (*_delta3_index)[e] = _delta3->PRE_HEAP; deba@326: } deba@326: for (int i = 0; i < _blossom_num; ++i) { kpeter@581: (*_delta2_index)[i] = _delta2->PRE_HEAP; kpeter@581: (*_delta4_index)[i] = _delta4->PRE_HEAP; deba@326: } alpar@877: deba@870: _unmatched = _node_num; deba@870: deba@867: _delta1->clear(); deba@867: _delta2->clear(); deba@867: _delta3->clear(); deba@867: _delta4->clear(); deba@867: _blossom_set->clear(); deba@867: _tree_set->clear(); deba@326: deba@326: int index = 0; deba@326: for (NodeIt n(_graph); n != INVALID; ++n) { deba@326: Value max = 0; deba@326: for (OutArcIt e(_graph, n); e != INVALID; ++e) { deba@326: if (_graph.target(e) == n) continue; deba@326: if ((dualScale * _weight[e]) / 2 > max) { deba@326: max = (dualScale * _weight[e]) / 2; deba@326: } deba@326: } kpeter@581: (*_node_index)[n] = index; deba@867: (*_node_data)[index].heap_index.clear(); deba@867: (*_node_data)[index].heap.clear(); deba@326: (*_node_data)[index].pot = max; deba@326: _delta1->push(n, max); deba@326: int blossom = deba@326: _blossom_set->insert(n, std::numeric_limits::max()); deba@326: deba@326: _tree_set->insert(blossom); deba@326: deba@326: (*_blossom_data)[blossom].status = EVEN; deba@326: (*_blossom_data)[blossom].pred = INVALID; deba@326: (*_blossom_data)[blossom].next = INVALID; deba@326: (*_blossom_data)[blossom].pot = 0; deba@326: (*_blossom_data)[blossom].offset = 0; deba@326: ++index; deba@326: } deba@326: for (EdgeIt e(_graph); e != INVALID; ++e) { deba@326: int si = (*_node_index)[_graph.u(e)]; deba@326: int ti = (*_node_index)[_graph.v(e)]; deba@326: if (_graph.u(e) != _graph.v(e)) { deba@326: _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - deba@326: dualScale * _weight[e]) / 2); deba@326: } deba@326: } deba@326: } deba@326: deba@870: /// \brief Initialize the algorithm with fractional matching deba@870: /// deba@870: /// This function initializes the algorithm with a fractional deba@870: /// matching. This initialization is also called jumpstart heuristic. deba@870: void fractionalInit() { deba@870: createStructures(); deba@876: deba@876: _blossom_node_list.clear(); deba@876: _blossom_potential.clear(); alpar@877: deba@870: if (_fractional == 0) { deba@870: _fractional = new FractionalMatching(_graph, _weight, false); deba@870: } deba@870: _fractional->run(); deba@870: deba@870: for (ArcIt e(_graph); e != INVALID; ++e) { deba@870: (*_node_heap_index)[e] = BinHeap::PRE_HEAP; deba@870: } deba@870: for (NodeIt n(_graph); n != INVALID; ++n) { deba@870: (*_delta1_index)[n] = _delta1->PRE_HEAP; deba@870: } deba@870: for (EdgeIt e(_graph); e != INVALID; ++e) { deba@870: (*_delta3_index)[e] = _delta3->PRE_HEAP; deba@870: } deba@870: for (int i = 0; i < _blossom_num; ++i) { deba@870: (*_delta2_index)[i] = _delta2->PRE_HEAP; deba@870: (*_delta4_index)[i] = _delta4->PRE_HEAP; deba@870: } deba@870: deba@870: _unmatched = 0; deba@870: deba@876: _delta1->clear(); deba@876: _delta2->clear(); deba@876: _delta3->clear(); deba@876: _delta4->clear(); deba@876: _blossom_set->clear(); deba@876: _tree_set->clear(); deba@876: deba@870: int index = 0; deba@870: for (NodeIt n(_graph); n != INVALID; ++n) { deba@870: Value pot = _fractional->nodeValue(n); deba@870: (*_node_index)[n] = index; deba@870: (*_node_data)[index].pot = pot; deba@876: (*_node_data)[index].heap_index.clear(); deba@876: (*_node_data)[index].heap.clear(); deba@870: int blossom = deba@870: _blossom_set->insert(n, std::numeric_limits::max()); deba@870: deba@870: (*_blossom_data)[blossom].status = MATCHED; deba@870: (*_blossom_data)[blossom].pred = INVALID; deba@870: (*_blossom_data)[blossom].next = _fractional->matching(n); deba@870: if (_fractional->matching(n) == INVALID) { deba@870: (*_blossom_data)[blossom].base = n; deba@870: } deba@870: (*_blossom_data)[blossom].pot = 0; deba@870: (*_blossom_data)[blossom].offset = 0; deba@870: ++index; deba@870: } deba@870: deba@870: typename Graph::template NodeMap processed(_graph, false); deba@870: for (NodeIt n(_graph); n != INVALID; ++n) { deba@870: if (processed[n]) continue; deba@870: processed[n] = true; deba@870: if (_fractional->matching(n) == INVALID) continue; deba@870: int num = 1; deba@870: Node v = _graph.target(_fractional->matching(n)); deba@870: while (n != v) { deba@870: processed[v] = true; deba@870: v = _graph.target(_fractional->matching(v)); deba@870: ++num; deba@870: } deba@870: deba@870: if (num % 2 == 1) { deba@870: std::vector subblossoms(num); deba@870: deba@870: subblossoms[--num] = _blossom_set->find(n); deba@870: _delta1->push(n, _fractional->nodeValue(n)); deba@870: v = _graph.target(_fractional->matching(n)); deba@870: while (n != v) { deba@870: subblossoms[--num] = _blossom_set->find(v); deba@870: _delta1->push(v, _fractional->nodeValue(v)); alpar@877: v = _graph.target(_fractional->matching(v)); deba@870: } alpar@877: alpar@877: int surface = deba@870: _blossom_set->join(subblossoms.begin(), subblossoms.end()); deba@870: (*_blossom_data)[surface].status = EVEN; deba@870: (*_blossom_data)[surface].pred = INVALID; deba@870: (*_blossom_data)[surface].next = INVALID; deba@870: (*_blossom_data)[surface].pot = 0; deba@870: (*_blossom_data)[surface].offset = 0; alpar@877: deba@870: _tree_set->insert(surface); deba@870: ++_unmatched; deba@870: } deba@870: } deba@870: deba@870: for (EdgeIt e(_graph); e != INVALID; ++e) { deba@870: int si = (*_node_index)[_graph.u(e)]; deba@870: int sb = _blossom_set->find(_graph.u(e)); deba@870: int ti = (*_node_index)[_graph.v(e)]; deba@870: int tb = _blossom_set->find(_graph.v(e)); deba@870: if ((*_blossom_data)[sb].status == EVEN && deba@870: (*_blossom_data)[tb].status == EVEN && sb != tb) { deba@870: _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - deba@870: dualScale * _weight[e]) / 2); deba@870: } deba@870: } deba@870: deba@870: for (NodeIt n(_graph); n != INVALID; ++n) { deba@870: int nb = _blossom_set->find(n); deba@870: if ((*_blossom_data)[nb].status != MATCHED) continue; deba@870: int ni = (*_node_index)[n]; deba@870: deba@870: for (OutArcIt e(_graph, n); e != INVALID; ++e) { deba@870: Node v = _graph.target(e); deba@870: int vb = _blossom_set->find(v); deba@870: int vi = (*_node_index)[v]; deba@870: deba@870: Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - deba@870: dualScale * _weight[e]; deba@870: deba@870: if ((*_blossom_data)[vb].status == EVEN) { deba@870: deba@870: int vt = _tree_set->find(vb); deba@870: deba@870: typename std::map::iterator it = deba@870: (*_node_data)[ni].heap_index.find(vt); deba@870: deba@870: if (it != (*_node_data)[ni].heap_index.end()) { deba@870: if ((*_node_data)[ni].heap[it->second] > rw) { deba@870: (*_node_data)[ni].heap.replace(it->second, e); deba@870: (*_node_data)[ni].heap.decrease(e, rw); deba@870: it->second = e; deba@870: } deba@870: } else { deba@870: (*_node_data)[ni].heap.push(e, rw); deba@870: (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); deba@870: } deba@870: } deba@870: } alpar@877: deba@870: if (!(*_node_data)[ni].heap.empty()) { deba@870: _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); deba@870: _delta2->push(nb, _blossom_set->classPrio(nb)); deba@870: } deba@870: } deba@870: } deba@870: kpeter@590: /// \brief Start the algorithm deba@326: /// kpeter@590: /// This function starts the algorithm. kpeter@590: /// deba@870: /// \pre \ref init() or \ref fractionalInit() must be called deba@870: /// before using this function. deba@326: void start() { deba@326: enum OpType { deba@326: D1, D2, D3, D4 deba@326: }; deba@326: deba@870: while (_unmatched > 0) { deba@326: Value d1 = !_delta1->empty() ? deba@326: _delta1->prio() : std::numeric_limits::max(); deba@326: deba@326: Value d2 = !_delta2->empty() ? deba@326: _delta2->prio() : std::numeric_limits::max(); deba@326: deba@326: Value d3 = !_delta3->empty() ? deba@326: _delta3->prio() : std::numeric_limits::max(); deba@326: deba@326: Value d4 = !_delta4->empty() ? deba@326: _delta4->prio() : std::numeric_limits::max(); deba@326: deba@868: _delta_sum = d3; OpType ot = D3; deba@868: if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; } deba@326: if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } deba@326: if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } deba@326: deba@326: switch (ot) { deba@326: case D1: deba@326: { deba@326: Node n = _delta1->top(); deba@326: unmatchNode(n); deba@870: --_unmatched; deba@326: } deba@326: break; deba@326: case D2: deba@326: { deba@326: int blossom = _delta2->top(); deba@326: Node n = _blossom_set->classTop(blossom); deba@868: Arc a = (*_node_data)[(*_node_index)[n]].heap.top(); deba@868: if ((*_blossom_data)[blossom].next == INVALID) { deba@868: augmentOnArc(a); deba@870: --_unmatched; deba@868: } else { deba@868: extendOnArc(a); deba@868: } deba@326: } deba@326: break; deba@326: case D3: deba@326: { deba@326: Edge e = _delta3->top(); deba@326: deba@326: int left_blossom = _blossom_set->find(_graph.u(e)); deba@326: int right_blossom = _blossom_set->find(_graph.v(e)); deba@326: deba@326: if (left_blossom == right_blossom) { deba@326: _delta3->pop(); deba@326: } else { deba@868: int left_tree = _tree_set->find(left_blossom); deba@868: int right_tree = _tree_set->find(right_blossom); deba@326: deba@326: if (left_tree == right_tree) { deba@327: shrinkOnEdge(e, left_tree); deba@326: } else { deba@327: augmentOnEdge(e); deba@870: _unmatched -= 2; deba@326: } deba@326: } deba@326: } break; deba@326: case D4: deba@326: splitBlossom(_delta4->top()); deba@326: break; deba@326: } deba@326: } deba@326: extractMatching(); deba@326: } deba@326: kpeter@590: /// \brief Run the algorithm. deba@326: /// kpeter@590: /// This method runs the \c %MaxWeightedMatching algorithm. deba@326: /// deba@326: /// \note mwm.run() is just a shortcut of the following code. deba@326: /// \code deba@870: /// mwm.fractionalInit(); deba@326: /// mwm.start(); deba@326: /// \endcode deba@326: void run() { deba@870: fractionalInit(); deba@326: start(); deba@326: } deba@326: deba@326: /// @} deba@326: kpeter@590: /// \name Primal Solution deba@868: /// Functions to get the primal solution, i.e. the maximum weighted kpeter@590: /// matching.\n kpeter@590: /// Either \ref run() or \ref start() function should be called before kpeter@590: /// using them. deba@326: deba@326: /// @{ deba@326: kpeter@590: /// \brief Return the weight of the matching. deba@326: /// kpeter@590: /// This function returns the weight of the found matching. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. kpeter@593: Value matchingWeight() const { deba@326: Value sum = 0; deba@326: for (NodeIt n(_graph); n != INVALID; ++n) { deba@326: if ((*_matching)[n] != INVALID) { deba@326: sum += _weight[(*_matching)[n]]; deba@326: } deba@326: } deba@868: return sum / 2; deba@326: } deba@326: kpeter@590: /// \brief Return the size (cardinality) of the matching. deba@326: /// kpeter@590: /// This function returns the size (cardinality) of the found matching. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@327: int matchingSize() const { deba@327: int num = 0; deba@327: for (NodeIt n(_graph); n != INVALID; ++n) { deba@327: if ((*_matching)[n] != INVALID) { deba@327: ++num; deba@327: } deba@327: } deba@327: return num /= 2; deba@327: } deba@327: kpeter@590: /// \brief Return \c true if the given edge is in the matching. deba@327: /// deba@868: /// This function returns \c true if the given edge is in the found kpeter@590: /// matching. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@327: bool matching(const Edge& edge) const { deba@327: return edge == (*_matching)[_graph.u(edge)]; deba@326: } deba@326: kpeter@590: /// \brief Return the matching arc (or edge) incident to the given node. deba@326: /// kpeter@590: /// This function returns the matching arc (or edge) incident to the deba@868: /// given node in the found matching or \c INVALID if the node is kpeter@590: /// not covered by the matching. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@326: Arc matching(const Node& node) const { deba@326: return (*_matching)[node]; deba@326: } deba@326: kpeter@593: /// \brief Return a const reference to the matching map. kpeter@593: /// kpeter@593: /// This function returns a const reference to a node map that stores kpeter@593: /// the matching arc (or edge) incident to each node. kpeter@593: const MatchingMap& matchingMap() const { kpeter@593: return *_matching; kpeter@593: } kpeter@593: kpeter@590: /// \brief Return the mate of the given node. deba@326: /// deba@868: /// This function returns the mate of the given node in the found kpeter@590: /// matching or \c INVALID if the node is not covered by the matching. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@326: Node mate(const Node& node) const { deba@326: return (*_matching)[node] != INVALID ? deba@326: _graph.target((*_matching)[node]) : INVALID; deba@326: } deba@326: deba@326: /// @} deba@326: kpeter@590: /// \name Dual Solution kpeter@590: /// Functions to get the dual solution.\n kpeter@590: /// Either \ref run() or \ref start() function should be called before kpeter@590: /// using them. deba@326: deba@326: /// @{ deba@326: kpeter@590: /// \brief Return the value of the dual solution. deba@326: /// deba@868: /// This function returns the value of the dual solution. deba@868: /// It should be equal to the primal value scaled by \ref dualScale kpeter@590: /// "dual scale". kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@326: Value dualValue() const { deba@326: Value sum = 0; deba@326: for (NodeIt n(_graph); n != INVALID; ++n) { deba@326: sum += nodeValue(n); deba@326: } deba@326: for (int i = 0; i < blossomNum(); ++i) { deba@326: sum += blossomValue(i) * (blossomSize(i) / 2); deba@326: } deba@326: return sum; deba@326: } deba@326: kpeter@590: /// \brief Return the dual value (potential) of the given node. deba@326: /// kpeter@590: /// This function returns the dual value (potential) of the given node. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@326: Value nodeValue(const Node& n) const { deba@326: return (*_node_potential)[n]; deba@326: } deba@326: kpeter@590: /// \brief Return the number of the blossoms in the basis. deba@326: /// kpeter@590: /// This function returns the number of the blossoms in the basis. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@326: /// \see BlossomIt deba@326: int blossomNum() const { deba@326: return _blossom_potential.size(); deba@326: } deba@326: kpeter@590: /// \brief Return the number of the nodes in the given blossom. deba@326: /// kpeter@590: /// This function returns the number of the nodes in the given blossom. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. kpeter@590: /// \see BlossomIt deba@326: int blossomSize(int k) const { deba@326: return _blossom_potential[k].end - _blossom_potential[k].begin; deba@326: } deba@326: kpeter@590: /// \brief Return the dual value (ptential) of the given blossom. deba@326: /// kpeter@590: /// This function returns the dual value (ptential) of the given blossom. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@326: Value blossomValue(int k) const { deba@326: return _blossom_potential[k].value; deba@326: } deba@326: kpeter@590: /// \brief Iterator for obtaining the nodes of a blossom. deba@326: /// deba@868: /// This class provides an iterator for obtaining the nodes of the kpeter@590: /// given blossom. It lists a subset of the nodes. deba@868: /// Before using this iterator, you must allocate a kpeter@590: /// MaxWeightedMatching class and execute it. deba@326: class BlossomIt { deba@326: public: deba@326: deba@326: /// \brief Constructor. deba@326: /// kpeter@590: /// Constructor to get the nodes of the given variable. kpeter@590: /// deba@868: /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or deba@868: /// \ref MaxWeightedMatching::start() "algorithm.start()" must be kpeter@590: /// called before initializing this iterator. deba@326: BlossomIt(const MaxWeightedMatching& algorithm, int variable) deba@326: : _algorithm(&algorithm) deba@326: { deba@326: _index = _algorithm->_blossom_potential[variable].begin; deba@326: _last = _algorithm->_blossom_potential[variable].end; deba@326: } deba@326: kpeter@590: /// \brief Conversion to \c Node. deba@326: /// kpeter@590: /// Conversion to \c Node. deba@326: operator Node() const { deba@327: return _algorithm->_blossom_node_list[_index]; deba@326: } deba@326: deba@326: /// \brief Increment operator. deba@326: /// deba@326: /// Increment operator. deba@326: BlossomIt& operator++() { deba@326: ++_index; deba@326: return *this; deba@326: } deba@326: deba@327: /// \brief Validity checking deba@327: /// deba@327: /// Checks whether the iterator is invalid. deba@327: bool operator==(Invalid) const { return _index == _last; } deba@327: deba@327: /// \brief Validity checking deba@327: /// deba@327: /// Checks whether the iterator is valid. deba@327: bool operator!=(Invalid) const { return _index != _last; } deba@326: deba@326: private: deba@326: const MaxWeightedMatching* _algorithm; deba@326: int _last; deba@326: int _index; deba@326: }; deba@326: deba@326: /// @} deba@326: deba@326: }; deba@326: deba@326: /// \ingroup matching deba@326: /// deba@326: /// \brief Weighted perfect matching in general graphs deba@326: /// deba@326: /// This class provides an efficient implementation of Edmond's deba@327: /// maximum weighted perfect matching algorithm. The implementation deba@326: /// is based on extensive use of priority queues and provides kpeter@559: /// \f$O(nm\log n)\f$ time complexity. deba@326: /// deba@868: /// The maximum weighted perfect matching problem is to find a subset of deba@868: /// the edges in an undirected graph with maximum overall weight for which kpeter@590: /// each node has exactly one incident edge. kpeter@590: /// It can be formulated with the following linear program. deba@326: /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f] deba@327: /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} deba@327: \quad \forall B\in\mathcal{O}\f] */ deba@326: /// \f[x_e \ge 0\quad \forall e\in E\f] deba@326: /// \f[\max \sum_{e\in E}x_ew_e\f] deba@327: /// where \f$\delta(X)\f$ is the set of edges incident to a node in deba@327: /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in deba@327: /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality deba@327: /// subsets of the nodes. deba@326: /// deba@326: /// The algorithm calculates an optimal matching and a proof of the deba@326: /// optimality. The solution of the dual problem can be used to check deba@327: /// the result of the algorithm. The dual linear problem is the kpeter@590: /// following. deba@327: /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge deba@327: w_{uv} \quad \forall uv\in E\f] */ deba@326: /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] deba@327: /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} deba@327: \frac{\vert B \vert - 1}{2}z_B\f] */ deba@326: /// deba@868: /// The algorithm can be executed with the run() function. kpeter@590: /// After it the matching (the primal solution) and the dual solution deba@868: /// can be obtained using the query functions and the deba@868: /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, deba@868: /// which is able to iterate on the nodes of a blossom. kpeter@590: /// If the value type is integer, then the dual solution is multiplied kpeter@590: /// by \ref MaxWeightedMatching::dualScale "4". kpeter@590: /// kpeter@593: /// \tparam GR The undirected graph type the algorithm runs on. deba@868: /// \tparam WM The type edge weight map. The default type is kpeter@590: /// \ref concepts::Graph::EdgeMap "GR::EdgeMap". kpeter@590: #ifdef DOXYGEN kpeter@590: template kpeter@590: #else kpeter@559: template > kpeter@590: #endif deba@326: class MaxWeightedPerfectMatching { deba@326: public: deba@326: kpeter@590: /// The graph type of the algorithm kpeter@559: typedef GR Graph; kpeter@590: /// The type of the edge weight map kpeter@559: typedef WM WeightMap; kpeter@590: /// The value type of the edge weights deba@326: typedef typename WeightMap::Value Value; deba@326: deba@326: /// \brief Scaling factor for dual solution deba@326: /// deba@326: /// Scaling factor for dual solution, it is equal to 4 or 1 deba@326: /// according to the value type. deba@326: static const int dualScale = deba@326: std::numeric_limits::is_integer ? 4 : 1; deba@326: kpeter@593: /// The type of the matching map deba@326: typedef typename Graph::template NodeMap deba@326: MatchingMap; deba@326: deba@326: private: deba@326: deba@326: TEMPLATE_GRAPH_TYPEDEFS(Graph); deba@326: deba@326: typedef typename Graph::template NodeMap NodePotential; deba@326: typedef std::vector BlossomNodeList; deba@326: deba@326: struct BlossomVariable { deba@326: int begin, end; deba@326: Value value; deba@326: deba@326: BlossomVariable(int _begin, int _end, Value _value) deba@326: : begin(_begin), end(_end), value(_value) {} deba@326: deba@326: }; deba@326: deba@326: typedef std::vector BlossomPotential; deba@326: deba@326: const Graph& _graph; deba@326: const WeightMap& _weight; deba@326: deba@326: MatchingMap* _matching; deba@326: deba@326: NodePotential* _node_potential; deba@326: deba@326: BlossomPotential _blossom_potential; deba@326: BlossomNodeList _blossom_node_list; deba@326: deba@326: int _node_num; deba@326: int _blossom_num; deba@326: deba@326: typedef RangeMap IntIntMap; deba@326: deba@326: enum Status { deba@326: EVEN = -1, MATCHED = 0, ODD = 1 deba@326: }; deba@326: deba@327: typedef HeapUnionFind BlossomSet; deba@326: struct BlossomData { deba@326: int tree; deba@326: Status status; deba@326: Arc pred, next; deba@326: Value pot, offset; deba@326: }; deba@326: deba@327: IntNodeMap *_blossom_index; deba@326: BlossomSet *_blossom_set; deba@326: RangeMap* _blossom_data; deba@326: deba@327: IntNodeMap *_node_index; deba@327: IntArcMap *_node_heap_index; deba@326: deba@326: struct NodeData { deba@326: deba@327: NodeData(IntArcMap& node_heap_index) deba@326: : heap(node_heap_index) {} deba@326: deba@326: int blossom; deba@326: Value pot; deba@327: BinHeap heap; deba@326: std::map heap_index; deba@326: deba@326: int tree; deba@326: }; deba@326: deba@326: RangeMap* _node_data; deba@326: deba@326: typedef ExtendFindEnum TreeSet; deba@326: deba@326: IntIntMap *_tree_set_index; deba@326: TreeSet *_tree_set; deba@326: deba@326: IntIntMap *_delta2_index; deba@326: BinHeap *_delta2; deba@326: deba@327: IntEdgeMap *_delta3_index; deba@327: BinHeap *_delta3; deba@326: deba@326: IntIntMap *_delta4_index; deba@326: BinHeap *_delta4; deba@326: deba@326: Value _delta_sum; deba@870: int _unmatched; deba@870: alpar@877: typedef MaxWeightedPerfectFractionalMatching deba@870: FractionalMatching; deba@870: FractionalMatching *_fractional; deba@326: deba@326: void createStructures() { deba@326: _node_num = countNodes(_graph); deba@326: _blossom_num = _node_num * 3 / 2; deba@326: deba@326: if (!_matching) { deba@326: _matching = new MatchingMap(_graph); deba@326: } deba@867: deba@326: if (!_node_potential) { deba@326: _node_potential = new NodePotential(_graph); deba@326: } deba@867: deba@326: if (!_blossom_set) { deba@327: _blossom_index = new IntNodeMap(_graph); deba@326: _blossom_set = new BlossomSet(*_blossom_index); deba@326: _blossom_data = new RangeMap(_blossom_num); deba@867: } else if (_blossom_data->size() != _blossom_num) { deba@867: delete _blossom_data; deba@867: _blossom_data = new RangeMap(_blossom_num); deba@326: } deba@326: deba@326: if (!_node_index) { deba@327: _node_index = new IntNodeMap(_graph); deba@327: _node_heap_index = new IntArcMap(_graph); deba@326: _node_data = new RangeMap(_node_num, deba@327: NodeData(*_node_heap_index)); deba@867: } else if (_node_data->size() != _node_num) { deba@867: delete _node_data; deba@867: _node_data = new RangeMap(_node_num, deba@867: NodeData(*_node_heap_index)); deba@326: } deba@326: deba@326: if (!_tree_set) { deba@326: _tree_set_index = new IntIntMap(_blossom_num); deba@326: _tree_set = new TreeSet(*_tree_set_index); deba@867: } else { deba@867: _tree_set_index->resize(_blossom_num); deba@326: } deba@867: deba@326: if (!_delta2) { deba@326: _delta2_index = new IntIntMap(_blossom_num); deba@326: _delta2 = new BinHeap(*_delta2_index); deba@867: } else { deba@867: _delta2_index->resize(_blossom_num); deba@326: } deba@867: deba@326: if (!_delta3) { deba@327: _delta3_index = new IntEdgeMap(_graph); deba@327: _delta3 = new BinHeap(*_delta3_index); deba@326: } deba@867: deba@326: if (!_delta4) { deba@326: _delta4_index = new IntIntMap(_blossom_num); deba@326: _delta4 = new BinHeap(*_delta4_index); deba@867: } else { deba@867: _delta4_index->resize(_blossom_num); deba@326: } deba@326: } deba@326: deba@326: void destroyStructures() { deba@326: if (_matching) { deba@326: delete _matching; deba@326: } deba@326: if (_node_potential) { deba@326: delete _node_potential; deba@326: } deba@326: if (_blossom_set) { deba@326: delete _blossom_index; deba@326: delete _blossom_set; deba@326: delete _blossom_data; deba@326: } deba@326: deba@326: if (_node_index) { deba@326: delete _node_index; deba@326: delete _node_heap_index; deba@326: delete _node_data; deba@326: } deba@326: deba@326: if (_tree_set) { deba@326: delete _tree_set_index; deba@326: delete _tree_set; deba@326: } deba@326: if (_delta2) { deba@326: delete _delta2_index; deba@326: delete _delta2; deba@326: } deba@326: if (_delta3) { deba@326: delete _delta3_index; deba@326: delete _delta3; deba@326: } deba@326: if (_delta4) { deba@326: delete _delta4_index; deba@326: delete _delta4; deba@326: } deba@326: } deba@326: deba@326: void matchedToEven(int blossom, int tree) { deba@326: if (_delta2->state(blossom) == _delta2->IN_HEAP) { deba@326: _delta2->erase(blossom); deba@326: } deba@326: deba@326: if (!_blossom_set->trivial(blossom)) { deba@326: (*_blossom_data)[blossom].pot -= deba@326: 2 * (_delta_sum - (*_blossom_data)[blossom].offset); deba@326: } deba@326: deba@326: for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); deba@326: n != INVALID; ++n) { deba@326: deba@326: _blossom_set->increase(n, std::numeric_limits::max()); deba@326: int ni = (*_node_index)[n]; deba@326: deba@326: (*_node_data)[ni].heap.clear(); deba@326: (*_node_data)[ni].heap_index.clear(); deba@326: deba@326: (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset; deba@326: deba@326: for (InArcIt e(_graph, n); e != INVALID; ++e) { deba@326: Node v = _graph.source(e); deba@326: int vb = _blossom_set->find(v); deba@326: int vi = (*_node_index)[v]; deba@326: deba@326: Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - deba@326: dualScale * _weight[e]; deba@326: deba@326: if ((*_blossom_data)[vb].status == EVEN) { deba@326: if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { deba@326: _delta3->push(e, rw / 2); deba@326: } deba@326: } else { deba@326: typename std::map::iterator it = deba@326: (*_node_data)[vi].heap_index.find(tree); deba@326: deba@326: if (it != (*_node_data)[vi].heap_index.end()) { deba@326: if ((*_node_data)[vi].heap[it->second] > rw) { deba@326: (*_node_data)[vi].heap.replace(it->second, e); deba@326: (*_node_data)[vi].heap.decrease(e, rw); deba@326: it->second = e; deba@326: } deba@326: } else { deba@326: (*_node_data)[vi].heap.push(e, rw); deba@326: (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); deba@326: } deba@326: deba@326: if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { deba@326: _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); deba@326: deba@326: if ((*_blossom_data)[vb].status == MATCHED) { deba@326: if (_delta2->state(vb) != _delta2->IN_HEAP) { deba@326: _delta2->push(vb, _blossom_set->classPrio(vb) - deba@326: (*_blossom_data)[vb].offset); deba@326: } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - deba@326: (*_blossom_data)[vb].offset){ deba@326: _delta2->decrease(vb, _blossom_set->classPrio(vb) - deba@326: (*_blossom_data)[vb].offset); deba@326: } deba@326: } deba@326: } deba@326: } deba@326: } deba@326: } deba@326: (*_blossom_data)[blossom].offset = 0; deba@326: } deba@326: deba@326: void matchedToOdd(int blossom) { deba@326: if (_delta2->state(blossom) == _delta2->IN_HEAP) { deba@326: _delta2->erase(blossom); deba@326: } deba@326: (*_blossom_data)[blossom].offset += _delta_sum; deba@326: if (!_blossom_set->trivial(blossom)) { deba@326: _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + deba@326: (*_blossom_data)[blossom].offset); deba@326: } deba@326: } deba@326: deba@326: void evenToMatched(int blossom, int tree) { deba@326: if (!_blossom_set->trivial(blossom)) { deba@326: (*_blossom_data)[blossom].pot += 2 * _delta_sum; deba@326: } deba@326: deba@326: for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); deba@326: n != INVALID; ++n) { deba@326: int ni = (*_node_index)[n]; deba@326: (*_node_data)[ni].pot -= _delta_sum; deba@326: deba@326: for (InArcIt e(_graph, n); e != INVALID; ++e) { deba@326: Node v = _graph.source(e); deba@326: int vb = _blossom_set->find(v); deba@326: int vi = (*_node_index)[v]; deba@326: deba@326: Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - deba@326: dualScale * _weight[e]; deba@326: deba@326: if (vb == blossom) { deba@326: if (_delta3->state(e) == _delta3->IN_HEAP) { deba@326: _delta3->erase(e); deba@326: } deba@326: } else if ((*_blossom_data)[vb].status == EVEN) { deba@326: deba@326: if (_delta3->state(e) == _delta3->IN_HEAP) { deba@326: _delta3->erase(e); deba@326: } deba@326: deba@326: int vt = _tree_set->find(vb); deba@326: deba@326: if (vt != tree) { deba@326: deba@326: Arc r = _graph.oppositeArc(e); deba@326: deba@326: typename std::map::iterator it = deba@326: (*_node_data)[ni].heap_index.find(vt); deba@326: deba@326: if (it != (*_node_data)[ni].heap_index.end()) { deba@326: if ((*_node_data)[ni].heap[it->second] > rw) { deba@326: (*_node_data)[ni].heap.replace(it->second, r); deba@326: (*_node_data)[ni].heap.decrease(r, rw); deba@326: it->second = r; deba@326: } deba@326: } else { deba@326: (*_node_data)[ni].heap.push(r, rw); deba@326: (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); deba@326: } deba@326: deba@326: if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { deba@326: _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); deba@326: deba@326: if (_delta2->state(blossom) != _delta2->IN_HEAP) { deba@326: _delta2->push(blossom, _blossom_set->classPrio(blossom) - deba@326: (*_blossom_data)[blossom].offset); deba@326: } else if ((*_delta2)[blossom] > deba@326: _blossom_set->classPrio(blossom) - deba@326: (*_blossom_data)[blossom].offset){ deba@326: _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - deba@326: (*_blossom_data)[blossom].offset); deba@326: } deba@326: } deba@326: } deba@326: } else { deba@326: deba@326: typename std::map::iterator it = deba@326: (*_node_data)[vi].heap_index.find(tree); deba@326: deba@326: if (it != (*_node_data)[vi].heap_index.end()) { deba@326: (*_node_data)[vi].heap.erase(it->second); deba@326: (*_node_data)[vi].heap_index.erase(it); deba@326: if ((*_node_data)[vi].heap.empty()) { deba@326: _blossom_set->increase(v, std::numeric_limits::max()); deba@326: } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { deba@326: _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); deba@326: } deba@326: deba@326: if ((*_blossom_data)[vb].status == MATCHED) { deba@326: if (_blossom_set->classPrio(vb) == deba@326: std::numeric_limits::max()) { deba@326: _delta2->erase(vb); deba@326: } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - deba@326: (*_blossom_data)[vb].offset) { deba@326: _delta2->increase(vb, _blossom_set->classPrio(vb) - deba@326: (*_blossom_data)[vb].offset); deba@326: } deba@326: } deba@326: } deba@326: } deba@326: } deba@326: } deba@326: } deba@326: deba@326: void oddToMatched(int blossom) { deba@326: (*_blossom_data)[blossom].offset -= _delta_sum; deba@326: deba@326: if (_blossom_set->classPrio(blossom) != deba@326: std::numeric_limits::max()) { deba@326: _delta2->push(blossom, _blossom_set->classPrio(blossom) - deba@326: (*_blossom_data)[blossom].offset); deba@326: } deba@326: deba@326: if (!_blossom_set->trivial(blossom)) { deba@326: _delta4->erase(blossom); deba@326: } deba@326: } deba@326: deba@326: void oddToEven(int blossom, int tree) { deba@326: if (!_blossom_set->trivial(blossom)) { deba@326: _delta4->erase(blossom); deba@326: (*_blossom_data)[blossom].pot -= deba@326: 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); deba@326: } deba@326: deba@326: for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); deba@326: n != INVALID; ++n) { deba@326: int ni = (*_node_index)[n]; deba@326: deba@326: _blossom_set->increase(n, std::numeric_limits::max()); deba@326: deba@326: (*_node_data)[ni].heap.clear(); deba@326: (*_node_data)[ni].heap_index.clear(); deba@326: (*_node_data)[ni].pot += deba@326: 2 * _delta_sum - (*_blossom_data)[blossom].offset; deba@326: deba@326: for (InArcIt e(_graph, n); e != INVALID; ++e) { deba@326: Node v = _graph.source(e); deba@326: int vb = _blossom_set->find(v); deba@326: int vi = (*_node_index)[v]; deba@326: deba@326: Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - deba@326: dualScale * _weight[e]; deba@326: deba@326: if ((*_blossom_data)[vb].status == EVEN) { deba@326: if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { deba@326: _delta3->push(e, rw / 2); deba@326: } deba@326: } else { deba@326: deba@326: typename std::map::iterator it = deba@326: (*_node_data)[vi].heap_index.find(tree); deba@326: deba@326: if (it != (*_node_data)[vi].heap_index.end()) { deba@326: if ((*_node_data)[vi].heap[it->second] > rw) { deba@326: (*_node_data)[vi].heap.replace(it->second, e); deba@326: (*_node_data)[vi].heap.decrease(e, rw); deba@326: it->second = e; deba@326: } deba@326: } else { deba@326: (*_node_data)[vi].heap.push(e, rw); deba@326: (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); deba@326: } deba@326: deba@326: if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { deba@326: _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); deba@326: deba@326: if ((*_blossom_data)[vb].status == MATCHED) { deba@326: if (_delta2->state(vb) != _delta2->IN_HEAP) { deba@326: _delta2->push(vb, _blossom_set->classPrio(vb) - deba@326: (*_blossom_data)[vb].offset); deba@326: } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - deba@326: (*_blossom_data)[vb].offset) { deba@326: _delta2->decrease(vb, _blossom_set->classPrio(vb) - deba@326: (*_blossom_data)[vb].offset); deba@326: } deba@326: } deba@326: } deba@326: } deba@326: } deba@326: } deba@326: (*_blossom_data)[blossom].offset = 0; deba@326: } deba@326: deba@326: void alternatePath(int even, int tree) { deba@326: int odd; deba@326: deba@326: evenToMatched(even, tree); deba@326: (*_blossom_data)[even].status = MATCHED; deba@326: deba@326: while ((*_blossom_data)[even].pred != INVALID) { deba@326: odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred)); deba@326: (*_blossom_data)[odd].status = MATCHED; deba@326: oddToMatched(odd); deba@326: (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; deba@326: deba@326: even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred)); deba@326: (*_blossom_data)[even].status = MATCHED; deba@326: evenToMatched(even, tree); deba@326: (*_blossom_data)[even].next = deba@326: _graph.oppositeArc((*_blossom_data)[odd].pred); deba@326: } deba@326: deba@326: } deba@326: deba@326: void destroyTree(int tree) { deba@326: for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { deba@326: if ((*_blossom_data)[b].status == EVEN) { deba@326: (*_blossom_data)[b].status = MATCHED; deba@326: evenToMatched(b, tree); deba@326: } else if ((*_blossom_data)[b].status == ODD) { deba@326: (*_blossom_data)[b].status = MATCHED; deba@326: oddToMatched(b); deba@326: } deba@326: } deba@326: _tree_set->eraseClass(tree); deba@326: } deba@326: deba@327: void augmentOnEdge(const Edge& edge) { deba@327: deba@327: int left = _blossom_set->find(_graph.u(edge)); deba@327: int right = _blossom_set->find(_graph.v(edge)); deba@326: deba@326: int left_tree = _tree_set->find(left); deba@326: alternatePath(left, left_tree); deba@326: destroyTree(left_tree); deba@326: deba@326: int right_tree = _tree_set->find(right); deba@326: alternatePath(right, right_tree); deba@326: destroyTree(right_tree); deba@326: deba@327: (*_blossom_data)[left].next = _graph.direct(edge, true); deba@327: (*_blossom_data)[right].next = _graph.direct(edge, false); deba@326: } deba@326: deba@326: void extendOnArc(const Arc& arc) { deba@326: int base = _blossom_set->find(_graph.target(arc)); deba@326: int tree = _tree_set->find(base); deba@326: deba@326: int odd = _blossom_set->find(_graph.source(arc)); deba@326: _tree_set->insert(odd, tree); deba@326: (*_blossom_data)[odd].status = ODD; deba@326: matchedToOdd(odd); deba@326: (*_blossom_data)[odd].pred = arc; deba@326: deba@326: int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next)); deba@326: (*_blossom_data)[even].pred = (*_blossom_data)[even].next; deba@326: _tree_set->insert(even, tree); deba@326: (*_blossom_data)[even].status = EVEN; deba@326: matchedToEven(even, tree); deba@326: } deba@326: deba@327: void shrinkOnEdge(const Edge& edge, int tree) { deba@326: int nca = -1; deba@326: std::vector left_path, right_path; deba@326: deba@326: { deba@326: std::set left_set, right_set; deba@326: int left = _blossom_set->find(_graph.u(edge)); deba@326: left_path.push_back(left); deba@326: left_set.insert(left); deba@326: deba@326: int right = _blossom_set->find(_graph.v(edge)); deba@326: right_path.push_back(right); deba@326: right_set.insert(right); deba@326: deba@326: while (true) { deba@326: deba@326: if ((*_blossom_data)[left].pred == INVALID) break; deba@326: deba@326: left = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); deba@326: left_path.push_back(left); deba@326: left = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); deba@326: left_path.push_back(left); deba@326: deba@326: left_set.insert(left); deba@326: deba@326: if (right_set.find(left) != right_set.end()) { deba@326: nca = left; deba@326: break; deba@326: } deba@326: deba@326: if ((*_blossom_data)[right].pred == INVALID) break; deba@326: deba@326: right = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); deba@326: right_path.push_back(right); deba@326: right = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); deba@326: right_path.push_back(right); deba@326: deba@326: right_set.insert(right); deba@326: deba@326: if (left_set.find(right) != left_set.end()) { deba@326: nca = right; deba@326: break; deba@326: } deba@326: deba@326: } deba@326: deba@326: if (nca == -1) { deba@326: if ((*_blossom_data)[left].pred == INVALID) { deba@326: nca = right; deba@326: while (left_set.find(nca) == left_set.end()) { deba@326: nca = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); deba@326: right_path.push_back(nca); deba@326: nca = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); deba@326: right_path.push_back(nca); deba@326: } deba@326: } else { deba@326: nca = left; deba@326: while (right_set.find(nca) == right_set.end()) { deba@326: nca = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); deba@326: left_path.push_back(nca); deba@326: nca = deba@326: _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); deba@326: left_path.push_back(nca); deba@326: } deba@326: } deba@326: } deba@326: } deba@326: deba@326: std::vector subblossoms; deba@326: Arc prev; deba@326: deba@326: prev = _graph.direct(edge, true); deba@326: for (int i = 0; left_path[i] != nca; i += 2) { deba@326: subblossoms.push_back(left_path[i]); deba@326: (*_blossom_data)[left_path[i]].next = prev; deba@326: _tree_set->erase(left_path[i]); deba@326: deba@326: subblossoms.push_back(left_path[i + 1]); deba@326: (*_blossom_data)[left_path[i + 1]].status = EVEN; deba@326: oddToEven(left_path[i + 1], tree); deba@326: _tree_set->erase(left_path[i + 1]); deba@326: prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred); deba@326: } deba@326: deba@326: int k = 0; deba@326: while (right_path[k] != nca) ++k; deba@326: deba@326: subblossoms.push_back(nca); deba@326: (*_blossom_data)[nca].next = prev; deba@326: deba@326: for (int i = k - 2; i >= 0; i -= 2) { deba@326: subblossoms.push_back(right_path[i + 1]); deba@326: (*_blossom_data)[right_path[i + 1]].status = EVEN; deba@326: oddToEven(right_path[i + 1], tree); deba@326: _tree_set->erase(right_path[i + 1]); deba@326: deba@326: (*_blossom_data)[right_path[i + 1]].next = deba@326: (*_blossom_data)[right_path[i + 1]].pred; deba@326: deba@326: subblossoms.push_back(right_path[i]); deba@326: _tree_set->erase(right_path[i]); deba@326: } deba@326: deba@326: int surface = deba@326: _blossom_set->join(subblossoms.begin(), subblossoms.end()); deba@326: deba@326: for (int i = 0; i < int(subblossoms.size()); ++i) { deba@326: if (!_blossom_set->trivial(subblossoms[i])) { deba@326: (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; deba@326: } deba@326: (*_blossom_data)[subblossoms[i]].status = MATCHED; deba@326: } deba@326: deba@326: (*_blossom_data)[surface].pot = -2 * _delta_sum; deba@326: (*_blossom_data)[surface].offset = 0; deba@326: (*_blossom_data)[surface].status = EVEN; deba@326: (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; deba@326: (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; deba@326: deba@326: _tree_set->insert(surface, tree); deba@326: _tree_set->erase(nca); deba@326: } deba@326: deba@326: void splitBlossom(int blossom) { deba@326: Arc next = (*_blossom_data)[blossom].next; deba@326: Arc pred = (*_blossom_data)[blossom].pred; deba@326: deba@326: int tree = _tree_set->find(blossom); deba@326: deba@326: (*_blossom_data)[blossom].status = MATCHED; deba@326: oddToMatched(blossom); deba@326: if (_delta2->state(blossom) == _delta2->IN_HEAP) { deba@326: _delta2->erase(blossom); deba@326: } deba@326: deba@326: std::vector subblossoms; deba@326: _blossom_set->split(blossom, std::back_inserter(subblossoms)); deba@326: deba@326: Value offset = (*_blossom_data)[blossom].offset; deba@326: int b = _blossom_set->find(_graph.source(pred)); deba@326: int d = _blossom_set->find(_graph.source(next)); deba@326: deba@326: int ib = -1, id = -1; deba@326: for (int i = 0; i < int(subblossoms.size()); ++i) { deba@326: if (subblossoms[i] == b) ib = i; deba@326: if (subblossoms[i] == d) id = i; deba@326: deba@326: (*_blossom_data)[subblossoms[i]].offset = offset; deba@326: if (!_blossom_set->trivial(subblossoms[i])) { deba@326: (*_blossom_data)[subblossoms[i]].pot -= 2 * offset; deba@326: } deba@326: if (_blossom_set->classPrio(subblossoms[i]) != deba@326: std::numeric_limits::max()) { deba@326: _delta2->push(subblossoms[i], deba@326: _blossom_set->classPrio(subblossoms[i]) - deba@326: (*_blossom_data)[subblossoms[i]].offset); deba@326: } deba@326: } deba@326: deba@326: if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) { deba@326: for (int i = (id + 1) % subblossoms.size(); deba@326: i != ib; i = (i + 2) % subblossoms.size()) { deba@326: int sb = subblossoms[i]; deba@326: int tb = subblossoms[(i + 1) % subblossoms.size()]; deba@326: (*_blossom_data)[sb].next = deba@326: _graph.oppositeArc((*_blossom_data)[tb].next); deba@326: } deba@326: deba@326: for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { deba@326: int sb = subblossoms[i]; deba@326: int tb = subblossoms[(i + 1) % subblossoms.size()]; deba@326: int ub = subblossoms[(i + 2) % subblossoms.size()]; deba@326: deba@326: (*_blossom_data)[sb].status = ODD; deba@326: matchedToOdd(sb); deba@326: _tree_set->insert(sb, tree); deba@326: (*_blossom_data)[sb].pred = pred; deba@326: (*_blossom_data)[sb].next = deba@326: _graph.oppositeArc((*_blossom_data)[tb].next); deba@326: deba@326: pred = (*_blossom_data)[ub].next; deba@326: deba@326: (*_blossom_data)[tb].status = EVEN; deba@326: matchedToEven(tb, tree); deba@326: _tree_set->insert(tb, tree); deba@326: (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; deba@326: } deba@326: deba@326: (*_blossom_data)[subblossoms[id]].status = ODD; deba@326: matchedToOdd(subblossoms[id]); deba@326: _tree_set->insert(subblossoms[id], tree); deba@326: (*_blossom_data)[subblossoms[id]].next = next; deba@326: (*_blossom_data)[subblossoms[id]].pred = pred; deba@326: deba@326: } else { deba@326: deba@326: for (int i = (ib + 1) % subblossoms.size(); deba@326: i != id; i = (i + 2) % subblossoms.size()) { deba@326: int sb = subblossoms[i]; deba@326: int tb = subblossoms[(i + 1) % subblossoms.size()]; deba@326: (*_blossom_data)[sb].next = deba@326: _graph.oppositeArc((*_blossom_data)[tb].next); deba@326: } deba@326: deba@326: for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { deba@326: int sb = subblossoms[i]; deba@326: int tb = subblossoms[(i + 1) % subblossoms.size()]; deba@326: int ub = subblossoms[(i + 2) % subblossoms.size()]; deba@326: deba@326: (*_blossom_data)[sb].status = ODD; deba@326: matchedToOdd(sb); deba@326: _tree_set->insert(sb, tree); deba@326: (*_blossom_data)[sb].next = next; deba@326: (*_blossom_data)[sb].pred = deba@326: _graph.oppositeArc((*_blossom_data)[tb].next); deba@326: deba@326: (*_blossom_data)[tb].status = EVEN; deba@326: matchedToEven(tb, tree); deba@326: _tree_set->insert(tb, tree); deba@326: (*_blossom_data)[tb].pred = deba@326: (*_blossom_data)[tb].next = deba@326: _graph.oppositeArc((*_blossom_data)[ub].next); deba@326: next = (*_blossom_data)[ub].next; deba@326: } deba@326: deba@326: (*_blossom_data)[subblossoms[ib]].status = ODD; deba@326: matchedToOdd(subblossoms[ib]); deba@326: _tree_set->insert(subblossoms[ib], tree); deba@326: (*_blossom_data)[subblossoms[ib]].next = next; deba@326: (*_blossom_data)[subblossoms[ib]].pred = pred; deba@326: } deba@326: _tree_set->erase(blossom); deba@326: } deba@326: deba@326: void extractBlossom(int blossom, const Node& base, const Arc& matching) { deba@326: if (_blossom_set->trivial(blossom)) { deba@326: int bi = (*_node_index)[base]; deba@326: Value pot = (*_node_data)[bi].pot; deba@326: kpeter@581: (*_matching)[base] = matching; deba@326: _blossom_node_list.push_back(base); kpeter@581: (*_node_potential)[base] = pot; deba@326: } else { deba@326: deba@326: Value pot = (*_blossom_data)[blossom].pot; deba@326: int bn = _blossom_node_list.size(); deba@326: deba@326: std::vector subblossoms; deba@326: _blossom_set->split(blossom, std::back_inserter(subblossoms)); deba@326: int b = _blossom_set->find(base); deba@326: int ib = -1; deba@326: for (int i = 0; i < int(subblossoms.size()); ++i) { deba@326: if (subblossoms[i] == b) { ib = i; break; } deba@326: } deba@326: deba@326: for (int i = 1; i < int(subblossoms.size()); i += 2) { deba@326: int sb = subblossoms[(ib + i) % subblossoms.size()]; deba@326: int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; deba@326: deba@326: Arc m = (*_blossom_data)[tb].next; deba@326: extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m)); deba@326: extractBlossom(tb, _graph.source(m), m); deba@326: } deba@326: extractBlossom(subblossoms[ib], base, matching); deba@326: deba@326: int en = _blossom_node_list.size(); deba@326: deba@326: _blossom_potential.push_back(BlossomVariable(bn, en, pot)); deba@326: } deba@326: } deba@326: deba@326: void extractMatching() { deba@326: std::vector blossoms; deba@326: for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { deba@326: blossoms.push_back(c); deba@326: } deba@326: deba@326: for (int i = 0; i < int(blossoms.size()); ++i) { deba@326: deba@326: Value offset = (*_blossom_data)[blossoms[i]].offset; deba@326: (*_blossom_data)[blossoms[i]].pot += 2 * offset; deba@326: for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); deba@326: n != INVALID; ++n) { deba@326: (*_node_data)[(*_node_index)[n]].pot -= offset; deba@326: } deba@326: deba@326: Arc matching = (*_blossom_data)[blossoms[i]].next; deba@326: Node base = _graph.source(matching); deba@326: extractBlossom(blossoms[i], base, matching); deba@326: } deba@326: } deba@326: deba@326: public: deba@326: deba@326: /// \brief Constructor deba@326: /// deba@326: /// Constructor. deba@326: MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight) deba@326: : _graph(graph), _weight(weight), _matching(0), deba@326: _node_potential(0), _blossom_potential(), _blossom_node_list(), deba@326: _node_num(0), _blossom_num(0), deba@326: deba@326: _blossom_index(0), _blossom_set(0), _blossom_data(0), deba@326: _node_index(0), _node_heap_index(0), _node_data(0), deba@326: _tree_set_index(0), _tree_set(0), deba@326: deba@326: _delta2_index(0), _delta2(0), deba@326: _delta3_index(0), _delta3(0), deba@326: _delta4_index(0), _delta4(0), deba@326: deba@870: _delta_sum(), _unmatched(0), deba@870: deba@870: _fractional(0) deba@870: {} deba@326: deba@326: ~MaxWeightedPerfectMatching() { deba@326: destroyStructures(); deba@870: if (_fractional) { deba@870: delete _fractional; deba@870: } deba@326: } deba@326: kpeter@590: /// \name Execution Control alpar@330: /// The simplest way to execute the algorithm is to use the kpeter@590: /// \ref run() member function. deba@326: deba@326: ///@{ deba@326: deba@326: /// \brief Initialize the algorithm deba@326: /// kpeter@590: /// This function initializes the algorithm. deba@326: void init() { deba@326: createStructures(); deba@326: deba@867: _blossom_node_list.clear(); deba@867: _blossom_potential.clear(); deba@867: deba@326: for (ArcIt e(_graph); e != INVALID; ++e) { kpeter@581: (*_node_heap_index)[e] = BinHeap::PRE_HEAP; deba@326: } deba@326: for (EdgeIt e(_graph); e != INVALID; ++e) { kpeter@581: (*_delta3_index)[e] = _delta3->PRE_HEAP; deba@326: } deba@326: for (int i = 0; i < _blossom_num; ++i) { kpeter@581: (*_delta2_index)[i] = _delta2->PRE_HEAP; kpeter@581: (*_delta4_index)[i] = _delta4->PRE_HEAP; deba@326: } deba@326: deba@870: _unmatched = _node_num; deba@870: deba@867: _delta2->clear(); deba@867: _delta3->clear(); deba@867: _delta4->clear(); deba@867: _blossom_set->clear(); deba@867: _tree_set->clear(); deba@867: deba@326: int index = 0; deba@326: for (NodeIt n(_graph); n != INVALID; ++n) { deba@326: Value max = - std::numeric_limits::max(); deba@326: for (OutArcIt e(_graph, n); e != INVALID; ++e) { deba@326: if (_graph.target(e) == n) continue; deba@326: if ((dualScale * _weight[e]) / 2 > max) { deba@326: max = (dualScale * _weight[e]) / 2; deba@326: } deba@326: } kpeter@581: (*_node_index)[n] = index; deba@867: (*_node_data)[index].heap_index.clear(); deba@867: (*_node_data)[index].heap.clear(); deba@326: (*_node_data)[index].pot = max; deba@326: int blossom = deba@326: _blossom_set->insert(n, std::numeric_limits::max()); deba@326: deba@326: _tree_set->insert(blossom); deba@326: deba@326: (*_blossom_data)[blossom].status = EVEN; deba@326: (*_blossom_data)[blossom].pred = INVALID; deba@326: (*_blossom_data)[blossom].next = INVALID; deba@326: (*_blossom_data)[blossom].pot = 0; deba@326: (*_blossom_data)[blossom].offset = 0; deba@326: ++index; deba@326: } deba@326: for (EdgeIt e(_graph); e != INVALID; ++e) { deba@326: int si = (*_node_index)[_graph.u(e)]; deba@326: int ti = (*_node_index)[_graph.v(e)]; deba@326: if (_graph.u(e) != _graph.v(e)) { deba@326: _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - deba@326: dualScale * _weight[e]) / 2); deba@326: } deba@326: } deba@326: } deba@326: deba@870: /// \brief Initialize the algorithm with fractional matching deba@870: /// deba@870: /// This function initializes the algorithm with a fractional deba@870: /// matching. This initialization is also called jumpstart heuristic. deba@870: void fractionalInit() { deba@870: createStructures(); deba@876: deba@876: _blossom_node_list.clear(); deba@876: _blossom_potential.clear(); alpar@877: deba@870: if (_fractional == 0) { deba@870: _fractional = new FractionalMatching(_graph, _weight, false); deba@870: } deba@870: if (!_fractional->run()) { deba@870: _unmatched = -1; deba@870: return; deba@870: } deba@870: deba@870: for (ArcIt e(_graph); e != INVALID; ++e) { deba@870: (*_node_heap_index)[e] = BinHeap::PRE_HEAP; deba@870: } deba@870: for (EdgeIt e(_graph); e != INVALID; ++e) { deba@870: (*_delta3_index)[e] = _delta3->PRE_HEAP; deba@870: } deba@870: for (int i = 0; i < _blossom_num; ++i) { deba@870: (*_delta2_index)[i] = _delta2->PRE_HEAP; deba@870: (*_delta4_index)[i] = _delta4->PRE_HEAP; deba@870: } deba@870: deba@870: _unmatched = 0; deba@870: deba@876: _delta2->clear(); deba@876: _delta3->clear(); deba@876: _delta4->clear(); deba@876: _blossom_set->clear(); deba@876: _tree_set->clear(); deba@876: deba@870: int index = 0; deba@870: for (NodeIt n(_graph); n != INVALID; ++n) { deba@870: Value pot = _fractional->nodeValue(n); deba@870: (*_node_index)[n] = index; deba@870: (*_node_data)[index].pot = pot; deba@876: (*_node_data)[index].heap_index.clear(); deba@876: (*_node_data)[index].heap.clear(); deba@870: int blossom = deba@870: _blossom_set->insert(n, std::numeric_limits::max()); deba@870: deba@870: (*_blossom_data)[blossom].status = MATCHED; deba@870: (*_blossom_data)[blossom].pred = INVALID; deba@870: (*_blossom_data)[blossom].next = _fractional->matching(n); deba@870: (*_blossom_data)[blossom].pot = 0; deba@870: (*_blossom_data)[blossom].offset = 0; deba@870: ++index; deba@870: } deba@870: deba@870: typename Graph::template NodeMap processed(_graph, false); deba@870: for (NodeIt n(_graph); n != INVALID; ++n) { deba@870: if (processed[n]) continue; deba@870: processed[n] = true; deba@870: if (_fractional->matching(n) == INVALID) continue; deba@870: int num = 1; deba@870: Node v = _graph.target(_fractional->matching(n)); deba@870: while (n != v) { deba@870: processed[v] = true; deba@870: v = _graph.target(_fractional->matching(v)); deba@870: ++num; deba@870: } deba@870: deba@870: if (num % 2 == 1) { deba@870: std::vector subblossoms(num); deba@870: deba@870: subblossoms[--num] = _blossom_set->find(n); deba@870: v = _graph.target(_fractional->matching(n)); deba@870: while (n != v) { deba@870: subblossoms[--num] = _blossom_set->find(v); alpar@877: v = _graph.target(_fractional->matching(v)); deba@870: } alpar@877: alpar@877: int surface = deba@870: _blossom_set->join(subblossoms.begin(), subblossoms.end()); deba@870: (*_blossom_data)[surface].status = EVEN; deba@870: (*_blossom_data)[surface].pred = INVALID; deba@870: (*_blossom_data)[surface].next = INVALID; deba@870: (*_blossom_data)[surface].pot = 0; deba@870: (*_blossom_data)[surface].offset = 0; alpar@877: deba@870: _tree_set->insert(surface); deba@870: ++_unmatched; deba@870: } deba@870: } deba@870: deba@870: for (EdgeIt e(_graph); e != INVALID; ++e) { deba@870: int si = (*_node_index)[_graph.u(e)]; deba@870: int sb = _blossom_set->find(_graph.u(e)); deba@870: int ti = (*_node_index)[_graph.v(e)]; deba@870: int tb = _blossom_set->find(_graph.v(e)); deba@870: if ((*_blossom_data)[sb].status == EVEN && deba@870: (*_blossom_data)[tb].status == EVEN && sb != tb) { deba@870: _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - deba@870: dualScale * _weight[e]) / 2); deba@870: } deba@870: } deba@870: deba@870: for (NodeIt n(_graph); n != INVALID; ++n) { deba@870: int nb = _blossom_set->find(n); deba@870: if ((*_blossom_data)[nb].status != MATCHED) continue; deba@870: int ni = (*_node_index)[n]; deba@870: deba@870: for (OutArcIt e(_graph, n); e != INVALID; ++e) { deba@870: Node v = _graph.target(e); deba@870: int vb = _blossom_set->find(v); deba@870: int vi = (*_node_index)[v]; deba@870: deba@870: Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - deba@870: dualScale * _weight[e]; deba@870: deba@870: if ((*_blossom_data)[vb].status == EVEN) { deba@870: deba@870: int vt = _tree_set->find(vb); deba@870: deba@870: typename std::map::iterator it = deba@870: (*_node_data)[ni].heap_index.find(vt); deba@870: deba@870: if (it != (*_node_data)[ni].heap_index.end()) { deba@870: if ((*_node_data)[ni].heap[it->second] > rw) { deba@870: (*_node_data)[ni].heap.replace(it->second, e); deba@870: (*_node_data)[ni].heap.decrease(e, rw); deba@870: it->second = e; deba@870: } deba@870: } else { deba@870: (*_node_data)[ni].heap.push(e, rw); deba@870: (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); deba@870: } deba@870: } deba@870: } alpar@877: deba@870: if (!(*_node_data)[ni].heap.empty()) { deba@870: _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); deba@870: _delta2->push(nb, _blossom_set->classPrio(nb)); deba@870: } deba@870: } deba@870: } deba@870: kpeter@590: /// \brief Start the algorithm deba@326: /// kpeter@590: /// This function starts the algorithm. kpeter@590: /// deba@870: /// \pre \ref init() or \ref fractionalInit() must be called before deba@870: /// using this function. deba@326: bool start() { deba@326: enum OpType { deba@326: D2, D3, D4 deba@326: }; deba@326: deba@870: if (_unmatched == -1) return false; deba@870: deba@870: while (_unmatched > 0) { deba@326: Value d2 = !_delta2->empty() ? deba@326: _delta2->prio() : std::numeric_limits::max(); deba@326: deba@326: Value d3 = !_delta3->empty() ? deba@326: _delta3->prio() : std::numeric_limits::max(); deba@326: deba@326: Value d4 = !_delta4->empty() ? deba@326: _delta4->prio() : std::numeric_limits::max(); deba@326: deba@868: _delta_sum = d3; OpType ot = D3; deba@868: if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } deba@326: if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } deba@326: deba@326: if (_delta_sum == std::numeric_limits::max()) { deba@326: return false; deba@326: } deba@326: deba@326: switch (ot) { deba@326: case D2: deba@326: { deba@326: int blossom = _delta2->top(); deba@326: Node n = _blossom_set->classTop(blossom); deba@326: Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); deba@326: extendOnArc(e); deba@326: } deba@326: break; deba@326: case D3: deba@326: { deba@326: Edge e = _delta3->top(); deba@326: deba@326: int left_blossom = _blossom_set->find(_graph.u(e)); deba@326: int right_blossom = _blossom_set->find(_graph.v(e)); deba@326: deba@326: if (left_blossom == right_blossom) { deba@326: _delta3->pop(); deba@326: } else { deba@326: int left_tree = _tree_set->find(left_blossom); deba@326: int right_tree = _tree_set->find(right_blossom); deba@326: deba@326: if (left_tree == right_tree) { deba@327: shrinkOnEdge(e, left_tree); deba@326: } else { deba@327: augmentOnEdge(e); deba@870: _unmatched -= 2; deba@326: } deba@326: } deba@326: } break; deba@326: case D4: deba@326: splitBlossom(_delta4->top()); deba@326: break; deba@326: } deba@326: } deba@326: extractMatching(); deba@326: return true; deba@326: } deba@326: kpeter@590: /// \brief Run the algorithm. deba@326: /// kpeter@590: /// This method runs the \c %MaxWeightedPerfectMatching algorithm. deba@326: /// kpeter@590: /// \note mwpm.run() is just a shortcut of the following code. deba@326: /// \code deba@870: /// mwpm.fractionalInit(); kpeter@590: /// mwpm.start(); deba@326: /// \endcode deba@326: bool run() { deba@870: fractionalInit(); deba@326: return start(); deba@326: } deba@326: deba@326: /// @} deba@326: kpeter@590: /// \name Primal Solution deba@868: /// Functions to get the primal solution, i.e. the maximum weighted kpeter@590: /// perfect matching.\n kpeter@590: /// Either \ref run() or \ref start() function should be called before kpeter@590: /// using them. deba@326: deba@326: /// @{ deba@326: kpeter@590: /// \brief Return the weight of the matching. deba@326: /// kpeter@590: /// This function returns the weight of the found matching. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. kpeter@593: Value matchingWeight() const { deba@326: Value sum = 0; deba@326: for (NodeIt n(_graph); n != INVALID; ++n) { deba@326: if ((*_matching)[n] != INVALID) { deba@326: sum += _weight[(*_matching)[n]]; deba@326: } deba@326: } deba@868: return sum / 2; deba@326: } deba@326: kpeter@590: /// \brief Return \c true if the given edge is in the matching. deba@326: /// deba@868: /// This function returns \c true if the given edge is in the found kpeter@590: /// matching. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@327: bool matching(const Edge& edge) const { deba@327: return static_cast((*_matching)[_graph.u(edge)]) == edge; deba@326: } deba@326: kpeter@590: /// \brief Return the matching arc (or edge) incident to the given node. deba@326: /// kpeter@590: /// This function returns the matching arc (or edge) incident to the deba@868: /// given node in the found matching or \c INVALID if the node is kpeter@590: /// not covered by the matching. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@326: Arc matching(const Node& node) const { deba@326: return (*_matching)[node]; deba@326: } deba@326: kpeter@593: /// \brief Return a const reference to the matching map. kpeter@593: /// kpeter@593: /// This function returns a const reference to a node map that stores kpeter@593: /// the matching arc (or edge) incident to each node. kpeter@593: const MatchingMap& matchingMap() const { kpeter@593: return *_matching; kpeter@593: } kpeter@593: kpeter@590: /// \brief Return the mate of the given node. deba@326: /// deba@868: /// This function returns the mate of the given node in the found kpeter@590: /// matching or \c INVALID if the node is not covered by the matching. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@326: Node mate(const Node& node) const { deba@326: return _graph.target((*_matching)[node]); deba@326: } deba@326: deba@326: /// @} deba@326: kpeter@590: /// \name Dual Solution kpeter@590: /// Functions to get the dual solution.\n kpeter@590: /// Either \ref run() or \ref start() function should be called before kpeter@590: /// using them. deba@326: deba@326: /// @{ deba@326: kpeter@590: /// \brief Return the value of the dual solution. deba@326: /// deba@868: /// This function returns the value of the dual solution. deba@868: /// It should be equal to the primal value scaled by \ref dualScale kpeter@590: /// "dual scale". kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@326: Value dualValue() const { deba@326: Value sum = 0; deba@326: for (NodeIt n(_graph); n != INVALID; ++n) { deba@326: sum += nodeValue(n); deba@326: } deba@326: for (int i = 0; i < blossomNum(); ++i) { deba@326: sum += blossomValue(i) * (blossomSize(i) / 2); deba@326: } deba@326: return sum; deba@326: } deba@326: kpeter@590: /// \brief Return the dual value (potential) of the given node. deba@326: /// kpeter@590: /// This function returns the dual value (potential) of the given node. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@326: Value nodeValue(const Node& n) const { deba@326: return (*_node_potential)[n]; deba@326: } deba@326: kpeter@590: /// \brief Return the number of the blossoms in the basis. deba@326: /// kpeter@590: /// This function returns the number of the blossoms in the basis. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@326: /// \see BlossomIt deba@326: int blossomNum() const { deba@326: return _blossom_potential.size(); deba@326: } deba@326: kpeter@590: /// \brief Return the number of the nodes in the given blossom. deba@326: /// kpeter@590: /// This function returns the number of the nodes in the given blossom. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. kpeter@590: /// \see BlossomIt deba@326: int blossomSize(int k) const { deba@326: return _blossom_potential[k].end - _blossom_potential[k].begin; deba@326: } deba@326: kpeter@590: /// \brief Return the dual value (ptential) of the given blossom. deba@326: /// kpeter@590: /// This function returns the dual value (ptential) of the given blossom. kpeter@590: /// kpeter@590: /// \pre Either run() or start() must be called before using this function. deba@326: Value blossomValue(int k) const { deba@326: return _blossom_potential[k].value; deba@326: } deba@326: kpeter@590: /// \brief Iterator for obtaining the nodes of a blossom. deba@326: /// deba@868: /// This class provides an iterator for obtaining the nodes of the kpeter@590: /// given blossom. It lists a subset of the nodes. deba@868: /// Before using this iterator, you must allocate a kpeter@590: /// MaxWeightedPerfectMatching class and execute it. deba@326: class BlossomIt { deba@326: public: deba@326: deba@326: /// \brief Constructor. deba@326: /// kpeter@590: /// Constructor to get the nodes of the given variable. kpeter@590: /// deba@868: /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" deba@868: /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" kpeter@590: /// must be called before initializing this iterator. deba@326: BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) deba@326: : _algorithm(&algorithm) deba@326: { deba@326: _index = _algorithm->_blossom_potential[variable].begin; deba@326: _last = _algorithm->_blossom_potential[variable].end; deba@326: } deba@326: kpeter@590: /// \brief Conversion to \c Node. deba@326: /// kpeter@590: /// Conversion to \c Node. deba@326: operator Node() const { deba@327: return _algorithm->_blossom_node_list[_index]; deba@326: } deba@326: deba@326: /// \brief Increment operator. deba@326: /// deba@326: /// Increment operator. deba@326: BlossomIt& operator++() { deba@326: ++_index; deba@326: return *this; deba@326: } deba@326: deba@327: /// \brief Validity checking deba@327: /// kpeter@590: /// This function checks whether the iterator is invalid. deba@327: bool operator==(Invalid) const { return _index == _last; } deba@327: deba@327: /// \brief Validity checking deba@327: /// kpeter@590: /// This function checks whether the iterator is valid. deba@327: bool operator!=(Invalid) const { return _index != _last; } deba@326: deba@326: private: deba@326: const MaxWeightedPerfectMatching* _algorithm; deba@326: int _last; deba@326: int _index; deba@326: }; deba@326: deba@326: /// @} deba@326: deba@326: }; deba@326: deba@326: } //END OF NAMESPACE LEMON deba@326: deba@868: #endif //LEMON_MATCHING_H