# HG changeset patch # User Balazs Dezso # Date 1223899211 -7200 # Node ID 91d63b8b1a4c6ba113680e8957b45d7610543417 # Parent 64ad48007fb28697c78f3f4c644894cd9572f14d Several improvements in maximum matching algorithms - The interface of MaxMatching is changed to be similar to the weighted algorithms - The internal data structure (the queue implementation and the matching map) is changed in the MaxMatching algorithm, which provides better runtime properties - The Blossom iterators are changed slightly in the weighted matching algorithms - Several documentation improvments - The test files are merged diff -r 64ad48007fb2 -r 91d63b8b1a4c lemon/max_matching.h --- a/lemon/max_matching.h Mon Oct 13 13:56:00 2008 +0200 +++ b/lemon/max_matching.h Mon Oct 13 14:00:11 2008 +0200 @@ -31,76 +31,405 @@ ///\ingroup matching ///\file -///\brief Maximum matching algorithms in graph. +///\brief Maximum matching algorithms in general graphs. namespace lemon { - ///\ingroup matching + /// \ingroup matching /// - ///\brief Edmonds' alternating forest maximum matching algorithm. + /// \brief Edmonds' alternating forest maximum matching algorithm. /// - ///This class provides Edmonds' alternating forest matching - ///algorithm. The starting matching (if any) can be passed to the - ///algorithm using some of init functions. + /// This class provides Edmonds' alternating forest matching + /// algorithm. The starting matching (if any) can be passed to the + /// algorithm using some of init functions. /// - ///The dual side of a matching is a map of the nodes to - ///MaxMatching::DecompType, having values \c D, \c A and \c C - ///showing the Gallai-Edmonds decomposition of the digraph. The nodes - ///in \c D induce a digraph with factor-critical components, the nodes - ///in \c A form the barrier, and the nodes in \c C induce a digraph - ///having a perfect matching. This decomposition can be attained by - ///calling \c decomposition() after running the algorithm. + /// The dual side of a matching is a map of the nodes to + /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c + /// MATCHED/C showing the Gallai-Edmonds decomposition of the + /// graph. The nodes in \c EVEN/D induce a graph with + /// factor-critical components, the nodes in \c ODD/A form the + /// barrier, and the nodes in \c MATCHED/C induce a graph having a + /// perfect matching. The number of the fractor critical components + /// minus the number of barrier nodes is a lower bound on the + /// unmatched nodes, and if the matching is optimal this bound is + /// tight. This decomposition can be attained by calling \c + /// decomposition() after running the algorithm. /// - ///\param Digraph The graph type the algorithm runs on. - template + /// \param _Graph The graph type the algorithm runs on. + template class MaxMatching { - - protected: + public: + + typedef _Graph Graph; + typedef typename Graph::template NodeMap + MatchingMap; + + ///\brief Indicates the Gallai-Edmonds decomposition of the graph. + /// + ///Indicates the Gallai-Edmonds decomposition of the graph, which + ///shows an upper bound on the size of a maximum matching. The + ///nodes with Status \c EVEN/D induce a graph with factor-critical + ///components, the nodes in \c ODD/A form the canonical barrier, + ///and the nodes in \c MATCHED/C induce a graph having a perfect + ///matching. + enum Status { + EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2 + }; + + typedef typename Graph::template NodeMap StatusMap; + + private: TEMPLATE_GRAPH_TYPEDEFS(Graph); - typedef typename Graph::template NodeMap UFECrossRef; - typedef UnionFindEnum UFE; - typedef std::vector NV; - - typedef typename Graph::template NodeMap EFECrossRef; - typedef ExtendFindEnum EFE; + typedef UnionFindEnum BlossomSet; + typedef ExtendFindEnum TreeSet; + typedef RangeMap NodeIntMap; + typedef MatchingMap EarMap; + typedef std::vector NodeQueue; + + const Graph& _graph; + MatchingMap* _matching; + StatusMap* _status; + + EarMap* _ear; + + IntNodeMap* _blossom_set_index; + BlossomSet* _blossom_set; + NodeIntMap* _blossom_rep; + + IntNodeMap* _tree_set_index; + TreeSet* _tree_set; + + NodeQueue _node_queue; + int _process, _postpone, _last; + + int _node_num; + + private: + + void createStructures() { + _node_num = countNodes(_graph); + if (!_matching) { + _matching = new MatchingMap(_graph); + } + if (!_status) { + _status = new StatusMap(_graph); + } + if (!_ear) { + _ear = new EarMap(_graph); + } + if (!_blossom_set) { + _blossom_set_index = new IntNodeMap(_graph); + _blossom_set = new BlossomSet(*_blossom_set_index); + } + if (!_blossom_rep) { + _blossom_rep = new NodeIntMap(_node_num); + } + if (!_tree_set) { + _tree_set_index = new IntNodeMap(_graph); + _tree_set = new TreeSet(*_tree_set_index); + } + _node_queue.resize(_node_num); + } + + void destroyStructures() { + if (_matching) { + delete _matching; + } + if (_status) { + delete _status; + } + if (_ear) { + delete _ear; + } + if (_blossom_set) { + delete _blossom_set; + delete _blossom_set_index; + } + if (_blossom_rep) { + delete _blossom_rep; + } + if (_tree_set) { + delete _tree_set_index; + delete _tree_set; + } + } + + void processDense(const Node& n) { + _process = _postpone = _last = 0; + _node_queue[_last++] = n; + + while (_process != _last) { + Node u = _node_queue[_process++]; + for (OutArcIt a(_graph, u); a != INVALID; ++a) { + Node v = _graph.target(a); + if ((*_status)[v] == MATCHED) { + extendOnArc(a); + } else if ((*_status)[v] == UNMATCHED) { + augmentOnArc(a); + return; + } + } + } + + while (_postpone != _last) { + Node u = _node_queue[_postpone++]; + + for (OutArcIt a(_graph, u); a != INVALID ; ++a) { + Node v = _graph.target(a); + + if ((*_status)[v] == EVEN) { + if (_blossom_set->find(u) != _blossom_set->find(v)) { + shrinkOnEdge(a); + } + } + + while (_process != _last) { + Node w = _node_queue[_process++]; + for (OutArcIt b(_graph, w); b != INVALID; ++b) { + Node x = _graph.target(b); + if ((*_status)[x] == MATCHED) { + extendOnArc(b); + } else if ((*_status)[x] == UNMATCHED) { + augmentOnArc(b); + return; + } + } + } + } + } + } + + void processSparse(const Node& n) { + _process = _last = 0; + _node_queue[_last++] = n; + while (_process != _last) { + Node u = _node_queue[_process++]; + for (OutArcIt a(_graph, u); a != INVALID; ++a) { + Node v = _graph.target(a); + + if ((*_status)[v] == EVEN) { + if (_blossom_set->find(u) != _blossom_set->find(v)) { + shrinkOnEdge(a); + } + } else if ((*_status)[v] == MATCHED) { + extendOnArc(a); + } else if ((*_status)[v] == UNMATCHED) { + augmentOnArc(a); + return; + } + } + } + } + + void shrinkOnEdge(const Edge& e) { + Node nca = INVALID; + + { + std::set left_set, right_set; + + Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))]; + left_set.insert(left); + + Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))]; + right_set.insert(right); + + while (true) { + if ((*_matching)[left] == INVALID) break; + left = _graph.target((*_matching)[left]); + left = (*_blossom_rep)[_blossom_set-> + find(_graph.target((*_ear)[left]))]; + if (right_set.find(left) != right_set.end()) { + nca = left; + break; + } + left_set.insert(left); + + if ((*_matching)[right] == INVALID) break; + right = _graph.target((*_matching)[right]); + right = (*_blossom_rep)[_blossom_set-> + find(_graph.target((*_ear)[right]))]; + if (left_set.find(right) != left_set.end()) { + nca = right; + break; + } + right_set.insert(right); + } + + if (nca == INVALID) { + if ((*_matching)[left] == INVALID) { + nca = right; + while (left_set.find(nca) == left_set.end()) { + nca = _graph.target((*_matching)[nca]); + nca =(*_blossom_rep)[_blossom_set-> + find(_graph.target((*_ear)[nca]))]; + } + } else { + nca = left; + while (right_set.find(nca) == right_set.end()) { + nca = _graph.target((*_matching)[nca]); + nca = (*_blossom_rep)[_blossom_set-> + find(_graph.target((*_ear)[nca]))]; + } + } + } + } + + { + + Node node = _graph.u(e); + Arc arc = _graph.direct(e, true); + Node base = (*_blossom_rep)[_blossom_set->find(node)]; + + while (base != nca) { + _ear->set(node, arc); + + Node n = node; + while (n != base) { + n = _graph.target((*_matching)[n]); + Arc a = (*_ear)[n]; + n = _graph.target(a); + _ear->set(n, _graph.oppositeArc(a)); + } + node = _graph.target((*_matching)[base]); + _tree_set->erase(base); + _tree_set->erase(node); + _blossom_set->insert(node, _blossom_set->find(base)); + _status->set(node, EVEN); + _node_queue[_last++] = node; + arc = _graph.oppositeArc((*_ear)[node]); + node = _graph.target((*_ear)[node]); + base = (*_blossom_rep)[_blossom_set->find(node)]; + _blossom_set->join(_graph.target(arc), base); + } + } + + _blossom_rep->set(_blossom_set->find(nca), nca); + + { + + Node node = _graph.v(e); + Arc arc = _graph.direct(e, false); + Node base = (*_blossom_rep)[_blossom_set->find(node)]; + + while (base != nca) { + _ear->set(node, arc); + + Node n = node; + while (n != base) { + n = _graph.target((*_matching)[n]); + Arc a = (*_ear)[n]; + n = _graph.target(a); + _ear->set(n, _graph.oppositeArc(a)); + } + node = _graph.target((*_matching)[base]); + _tree_set->erase(base); + _tree_set->erase(node); + _blossom_set->insert(node, _blossom_set->find(base)); + _status->set(node, EVEN); + _node_queue[_last++] = node; + arc = _graph.oppositeArc((*_ear)[node]); + node = _graph.target((*_ear)[node]); + base = (*_blossom_rep)[_blossom_set->find(node)]; + _blossom_set->join(_graph.target(arc), base); + } + } + + _blossom_rep->set(_blossom_set->find(nca), nca); + } + + + + void extendOnArc(const Arc& a) { + Node base = _graph.source(a); + Node odd = _graph.target(a); + + _ear->set(odd, _graph.oppositeArc(a)); + Node even = _graph.target((*_matching)[odd]); + _blossom_rep->set(_blossom_set->insert(even), even); + _status->set(odd, ODD); + _status->set(even, EVEN); + int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]); + _tree_set->insert(odd, tree); + _tree_set->insert(even, tree); + _node_queue[_last++] = even; + + } + + void augmentOnArc(const Arc& a) { + Node even = _graph.source(a); + Node odd = _graph.target(a); + + int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]); + + _matching->set(odd, _graph.oppositeArc(a)); + _status->set(odd, MATCHED); + + Arc arc = (*_matching)[even]; + _matching->set(even, a); + + while (arc != INVALID) { + odd = _graph.target(arc); + arc = (*_ear)[odd]; + even = _graph.target(arc); + _matching->set(odd, arc); + arc = (*_matching)[even]; + _matching->set(even, _graph.oppositeArc((*_matching)[odd])); + } + + for (typename TreeSet::ItemIt it(*_tree_set, tree); + it != INVALID; ++it) { + if ((*_status)[it] == ODD) { + _status->set(it, MATCHED); + } else { + int blossom = _blossom_set->find(it); + for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom); + jt != INVALID; ++jt) { + _status->set(jt, MATCHED); + } + _blossom_set->eraseClass(blossom); + } + } + _tree_set->eraseClass(tree); + + } public: - ///\brief Indicates the Gallai-Edmonds decomposition of the digraph. + /// \brief Constructor /// - ///Indicates the Gallai-Edmonds decomposition of the digraph, which - ///shows an upper bound on the size of a maximum matching. The - ///nodes with DecompType \c D induce a digraph with factor-critical - ///components, the nodes in \c A form the canonical barrier, and the - ///nodes in \c C induce a digraph having a perfect matching. - enum DecompType { - D=0, - A=1, - C=2 - }; - - protected: - - static const int HEUR_density=2; - const Graph& g; - typename Graph::template NodeMap _mate; - typename Graph::template NodeMap position; - - public: - - MaxMatching(const Graph& _g) - : g(_g), _mate(_g), position(_g) {} - - ///\brief Sets the actual matching to the empty matching. + /// Constructor. + MaxMatching(const Graph& graph) + : _graph(graph), _matching(0), _status(0), _ear(0), + _blossom_set_index(0), _blossom_set(0), _blossom_rep(0), + _tree_set_index(0), _tree_set(0) {} + + ~MaxMatching() { + destroyStructures(); + } + + /// \name Execution control + /// The simplest way to execute the algorithm is to use the member + /// \c run() member function. + /// \n + + /// If you need more control on the execution, first you must call + /// \ref init(), \ref greedyInit() or \ref matchingInit() + /// functions, then you can start the algorithm with the \ref + /// startParse() or startDense() functions. + + ///@{ + + /// \brief Sets the actual matching to the empty matching. /// - ///Sets the actual matching to the empty matching. + /// Sets the actual matching to the empty matching. /// void init() { - for(NodeIt v(g); v!=INVALID; ++v) { - _mate.set(v,INVALID); - position.set(v,C); + createStructures(); + for(NodeIt n(_graph); n != INVALID; ++n) { + _matching->set(n, INVALID); + _status->set(n, UNMATCHED); } } @@ -108,88 +437,96 @@ /// ///For initial matchig it finds a maximal greedy matching. void greedyInit() { - for(NodeIt v(g); v!=INVALID; ++v) { - _mate.set(v,INVALID); - position.set(v,C); + createStructures(); + for (NodeIt n(_graph); n != INVALID; ++n) { + _matching->set(n, INVALID); + _status->set(n, UNMATCHED); } - for(NodeIt v(g); v!=INVALID; ++v) - if ( _mate[v]==INVALID ) { - for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) { - Node y=g.runningNode(e); - if ( _mate[y]==INVALID && y!=v ) { - _mate.set(v,y); - _mate.set(y,v); + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_matching)[n] == INVALID) { + for (OutArcIt a(_graph, n); a != INVALID ; ++a) { + Node v = _graph.target(a); + if ((*_matching)[v] == INVALID && v != n) { + _matching->set(n, a); + _status->set(n, MATCHED); + _matching->set(v, _graph.oppositeArc(a)); + _status->set(v, MATCHED); break; } } } - } - - ///\brief Initialize the matching from each nodes' mate. - /// - ///Initialize the matching from a \c Node valued \c Node map. This - ///map must be \e symmetric, i.e. if \c map[u]==v then \c - ///map[v]==u must hold, and \c uv will be an arc of the initial - ///matching. - template - void mateMapInit(MateMap& map) { - for(NodeIt v(g); v!=INVALID; ++v) { - _mate.set(v,map[v]); - position.set(v,C); } } - ///\brief Initialize the matching from a node map with the - ///incident matching arcs. + + /// \brief Initialize the matching from the map containing a matching. /// - ///Initialize the matching from an \c Edge valued \c Node map. \c - ///map[v] must be an \c Edge incident to \c v. This map must have - ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c - ///g.oppositeNode(v,map[v])==u holds, and now some arc joining \c - ///u to \c v will be an arc of the matching. - template - void matchingMapInit(MatchingMap& map) { - for(NodeIt v(g); v!=INVALID; ++v) { - position.set(v,C); - Edge e=map[v]; - if ( e!=INVALID ) - _mate.set(v,g.oppositeNode(v,e)); - else - _mate.set(v,INVALID); + /// Initialize the matching from a \c bool valued \c Edge map. This + /// map must have the property that there are no two incident edges + /// with true value, ie. it contains a matching. + /// \return %True if the map contains a matching. + template + bool matchingInit(const MatchingMap& matching) { + createStructures(); + + for (NodeIt n(_graph); n != INVALID; ++n) { + _matching->set(n, INVALID); + _status->set(n, UNMATCHED); } + for(EdgeIt e(_graph); e!=INVALID; ++e) { + if (matching[e]) { + + Node u = _graph.u(e); + if ((*_matching)[u] != INVALID) return false; + _matching->set(u, _graph.direct(e, true)); + _status->set(u, MATCHED); + + Node v = _graph.v(e); + if ((*_matching)[v] != INVALID) return false; + _matching->set(v, _graph.direct(e, false)); + _status->set(v, MATCHED); + } + } + return true; } - ///\brief Initialize the matching from the map containing the - ///undirected matching arcs. + /// \brief Starts Edmonds' algorithm /// - ///Initialize the matching from a \c bool valued \c Edge map. This - ///map must have the property that there are no two incident arcs - ///\c e, \c f with \c map[e]==map[f]==true. The arcs \c e with \c - ///map[e]==true form the matching. - template - void matchingInit(MatchingMap& map) { - for(NodeIt v(g); v!=INVALID; ++v) { - _mate.set(v,INVALID); - position.set(v,C); - } - for(EdgeIt e(g); e!=INVALID; ++e) { - if ( map[e] ) { - Node u=g.u(e); - Node v=g.v(e); - _mate.set(u,v); - _mate.set(v,u); + /// If runs the original Edmonds' algorithm. + void startSparse() { + for(NodeIt n(_graph); n != INVALID; ++n) { + if ((*_status)[n] == UNMATCHED) { + (*_blossom_rep)[_blossom_set->insert(n)] = n; + _tree_set->insert(n); + _status->set(n, EVEN); + processSparse(n); } } } - - ///\brief Runs Edmonds' algorithm. + /// \brief Starts Edmonds' algorithm. /// - ///Runs Edmonds' algorithm for sparse digraphs (number of arcs < - ///2*number of nodes), and a heuristical Edmonds' algorithm with a - ///heuristic of postponing shrinks for dense digraphs. + /// It runs Edmonds' algorithm with a heuristic of postponing + /// shrinks, giving a faster algorithm for dense graphs. + void startDense() { + for(NodeIt n(_graph); n != INVALID; ++n) { + if ((*_status)[n] == UNMATCHED) { + (*_blossom_rep)[_blossom_set->insert(n)] = n; + _tree_set->insert(n); + _status->set(n, EVEN); + processDense(n); + } + } + } + + + /// \brief Runs Edmonds' algorithm + /// + /// Runs Edmonds' algorithm for sparse graphs (m<2*n) + /// or Edmonds' algorithm with a heuristic of + /// postponing shrinks for dense graphs. void run() { - if (countEdges(g) < HEUR_density * countNodes(g)) { + if (countEdges(_graph) < 2 * countNodes(_graph)) { greedyInit(); startSparse(); } else { @@ -198,422 +535,75 @@ } } - - ///\brief Starts Edmonds' algorithm. - /// - ///If runs the original Edmonds' algorithm. - void startSparse() { - - typename Graph::template NodeMap ear(g,INVALID); - //undefined for the base nodes of the blossoms (i.e. for the - //representative elements of UFE blossom) and for the nodes in C - - UFECrossRef blossom_base(g); - UFE blossom(blossom_base); - NV rep(countNodes(g)); - - EFECrossRef tree_base(g); - EFE tree(tree_base); - - //If these UFE's would be members of the class then also - //blossom_base and tree_base should be a member. - - //We build only one tree and the other vertices uncovered by the - //matching belong to C. (They can be considered as singleton - //trees.) If this tree can be augmented or no more - //grow/augmentation/shrink is possible then we return to this - //"for" cycle. - for(NodeIt v(g); v!=INVALID; ++v) { - if (position[v]==C && _mate[v]==INVALID) { - rep[blossom.insert(v)] = v; - tree.insert(v); - position.set(v,D); - normShrink(v, ear, blossom, rep, tree); - } - } - } - - ///\brief Starts Edmonds' algorithm. - /// - ///It runs Edmonds' algorithm with a heuristic of postponing - ///shrinks, giving a faster algorithm for dense digraphs. - void startDense() { - - typename Graph::template NodeMap ear(g,INVALID); - //undefined for the base nodes of the blossoms (i.e. for the - //representative elements of UFE blossom) and for the nodes in C - - UFECrossRef blossom_base(g); - UFE blossom(blossom_base); - NV rep(countNodes(g)); - - EFECrossRef tree_base(g); - EFE tree(tree_base); - - //If these UFE's would be members of the class then also - //blossom_base and tree_base should be a member. - - //We build only one tree and the other vertices uncovered by the - //matching belong to C. (They can be considered as singleton - //trees.) If this tree can be augmented or no more - //grow/augmentation/shrink is possible then we return to this - //"for" cycle. - for(NodeIt v(g); v!=INVALID; ++v) { - if ( position[v]==C && _mate[v]==INVALID ) { - rep[blossom.insert(v)] = v; - tree.insert(v); - position.set(v,D); - lateShrink(v, ear, blossom, rep, tree); - } - } - } - - + /// @} + + /// \name Primal solution + /// Functions for get the primal solution, ie. the matching. + + /// @{ ///\brief Returns the size of the actual matching stored. /// ///Returns the size of the actual matching stored. After \ref - ///run() it returns the size of a maximum matching in the digraph. - int size() const { - int s=0; - for(NodeIt v(g); v!=INVALID; ++v) { - if ( _mate[v]!=INVALID ) { - ++s; + ///run() it returns the size of the maximum matching in the graph. + int matchingSize() const { + int size = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_matching)[n] != INVALID) { + ++size; } } - return s/2; + return size / 2; } + /// \brief Returns true when the edge is in the matching. + /// + /// Returns true when the edge is in the matching. + bool matching(const Edge& edge) const { + return edge == (*_matching)[_graph.u(edge)]; + } + + /// \brief Returns the matching edge incident to the given node. + /// + /// Returns the matching edge of a \c node in the actual matching or + /// INVALID if the \c node is not covered by the actual matching. + Arc matching(const Node& n) const { + return (*_matching)[n]; + } ///\brief Returns the mate of a node in the actual matching. /// - ///Returns the mate of a \c node in the actual matching. - ///Returns INVALID if the \c node is not covered by the actual matching. - Node mate(const Node& node) const { - return _mate[node]; + ///Returns the mate of a \c node in the actual matching or + ///INVALID if the \c node is not covered by the actual matching. + Node mate(const Node& n) const { + return (*_matching)[n] != INVALID ? + _graph.target((*_matching)[n]) : INVALID; } - ///\brief Returns the matching arc incident to the given node. - /// - ///Returns the matching arc of a \c node in the actual matching. - ///Returns INVALID if the \c node is not covered by the actual matching. - Edge matchingArc(const Node& node) const { - if (_mate[node] == INVALID) return INVALID; - Node n = node < _mate[node] ? node : _mate[node]; - for (IncEdgeIt e(g, n); e != INVALID; ++e) { - if (g.oppositeNode(n, e) == _mate[n]) { - return e; - } - } - return INVALID; - } + /// @} + + /// \name Dual solution + /// Functions for get the dual solution, ie. the decomposition. + + /// @{ /// \brief Returns the class of the node in the Edmonds-Gallai /// decomposition. /// /// Returns the class of the node in the Edmonds-Gallai /// decomposition. - DecompType decomposition(const Node& n) { - return position[n] == A; + Status decomposition(const Node& n) const { + return (*_status)[n]; } /// \brief Returns true when the node is in the barrier. /// /// Returns true when the node is in the barrier. - bool barrier(const Node& n) { - return position[n] == A; + bool barrier(const Node& n) const { + return (*_status)[n] == ODD; } - ///\brief Gives back the matching in a \c Node of mates. - /// - ///Writes the stored matching to a \c Node valued \c Node map. The - ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c - ///map[v]==u will hold, and now \c uv is an arc of the matching. - template - void mateMap(MateMap& map) const { - for(NodeIt v(g); v!=INVALID; ++v) { - map.set(v,_mate[v]); - } - } - - ///\brief Gives back the matching in an \c Edge valued \c Node - ///map. - /// - ///Writes the stored matching to an \c Edge valued \c Node - ///map. \c map[v] will be an \c Edge incident to \c v. This - ///map will have the property that if \c g.oppositeNode(u,map[u]) - ///== v then \c map[u]==map[v] holds, and now this arc is an arc - ///of the matching. - template - void matchingMap(MatchingMap& map) const { - typename Graph::template NodeMap todo(g,true); - for(NodeIt v(g); v!=INVALID; ++v) { - if (_mate[v]!=INVALID && v < _mate[v]) { - Node u=_mate[v]; - for(IncEdgeIt e(g,v); e!=INVALID; ++e) { - if ( g.runningNode(e) == u ) { - map.set(u,e); - map.set(v,e); - todo.set(u,false); - todo.set(v,false); - break; - } - } - } - } - } - - - ///\brief Gives back the matching in a \c bool valued \c Edge - ///map. - /// - ///Writes the matching stored to a \c bool valued \c Arc - ///map. This map will have the property that there are no two - ///incident arcs \c e, \c f with \c map[e]==map[f]==true. The - ///arcs \c e with \c map[e]==true form the matching. - template - void matching(MatchingMap& map) const { - for(EdgeIt e(g); e!=INVALID; ++e) map.set(e,false); - - typename Graph::template NodeMap todo(g,true); - for(NodeIt v(g); v!=INVALID; ++v) { - if ( todo[v] && _mate[v]!=INVALID ) { - Node u=_mate[v]; - for(IncEdgeIt e(g,v); e!=INVALID; ++e) { - if ( g.runningNode(e) == u ) { - map.set(e,true); - todo.set(u,false); - todo.set(v,false); - break; - } - } - } - } - } - - - ///\brief Returns the canonical decomposition of the digraph after running - ///the algorithm. - /// - ///After calling any run methods of the class, it writes the - ///Gallai-Edmonds canonical decomposition of the digraph. \c map - ///must be a node map of \ref DecompType 's. - template - void decomposition(DecompositionMap& map) const { - for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]); - } - - ///\brief Returns a barrier on the nodes. - /// - ///After calling any run methods of the class, it writes a - ///canonical barrier on the nodes. The odd component number of the - ///remaining digraph minus the barrier size is a lower bound for the - ///uncovered nodes in the digraph. The \c map must be a node map of - ///bools. - template - void barrier(BarrierMap& barrier) { - for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A); - } - - private: - - - void lateShrink(Node v, typename Graph::template NodeMap& ear, - UFE& blossom, NV& rep, EFE& tree) { - //We have one tree which we grow, and also shrink but only if it - //cannot be postponed. If we augment then we return to the "for" - //cycle of runEdmonds(). - - std::queue Q; //queue of the totally unscanned nodes - Q.push(v); - std::queue R; - //queue of the nodes which must be scanned for a possible shrink - - while ( !Q.empty() ) { - Node x=Q.front(); - Q.pop(); - for( IncEdgeIt e(g,x); e!= INVALID; ++e ) { - Node y=g.runningNode(e); - //growOrAugment grows if y is covered by the matching and - //augments if not. In this latter case it returns 1. - if (position[y]==C && - growOrAugment(y, x, ear, blossom, rep, tree, Q)) return; - } - R.push(x); - } - - while ( !R.empty() ) { - Node x=R.front(); - R.pop(); - - for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) { - Node y=g.runningNode(e); - - if ( position[y] == D && blossom.find(x) != blossom.find(y) ) - //Recall that we have only one tree. - shrink( x, y, ear, blossom, rep, tree, Q); - - while ( !Q.empty() ) { - Node z=Q.front(); - Q.pop(); - for( IncEdgeIt f(g,z); f!= INVALID; ++f ) { - Node w=g.runningNode(f); - //growOrAugment grows if y is covered by the matching and - //augments if not. In this latter case it returns 1. - if (position[w]==C && - growOrAugment(w, z, ear, blossom, rep, tree, Q)) return; - } - R.push(z); - } - } //for e - } // while ( !R.empty() ) - } - - void normShrink(Node v, typename Graph::template NodeMap& ear, - UFE& blossom, NV& rep, EFE& tree) { - //We have one tree, which we grow and shrink. If we augment then we - //return to the "for" cycle of runEdmonds(). - - std::queue Q; //queue of the unscanned nodes - Q.push(v); - while ( !Q.empty() ) { - - Node x=Q.front(); - Q.pop(); - - for( IncEdgeIt e(g,x); e!=INVALID; ++e ) { - Node y=g.runningNode(e); - - switch ( position[y] ) { - case D: //x and y must be in the same tree - if ( blossom.find(x) != blossom.find(y)) - //x and y are in the same tree - shrink(x, y, ear, blossom, rep, tree, Q); - break; - case C: - //growOrAugment grows if y is covered by the matching and - //augments if not. In this latter case it returns 1. - if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return; - break; - default: break; - } - } - } - } - - void shrink(Node x,Node y, typename Graph::template NodeMap& ear, - UFE& blossom, NV& rep, EFE& tree,std::queue& Q) { - //x and y are the two adjacent vertices in two blossoms. - - typename Graph::template NodeMap path(g,false); - - Node b=rep[blossom.find(x)]; - path.set(b,true); - b=_mate[b]; - while ( b!=INVALID ) { - b=rep[blossom.find(ear[b])]; - path.set(b,true); - b=_mate[b]; - } //we go until the root through bases of blossoms and odd vertices - - Node top=y; - Node middle=rep[blossom.find(top)]; - Node bottom=x; - while ( !path[middle] ) - shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q); - //Until we arrive to a node on the path, we update blossom, tree - //and the positions of the odd nodes. - - Node base=middle; - top=x; - middle=rep[blossom.find(top)]; - bottom=y; - Node blossom_base=rep[blossom.find(base)]; - while ( middle!=blossom_base ) - shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q); - //Until we arrive to a node on the path, we update blossom, tree - //and the positions of the odd nodes. - - rep[blossom.find(base)] = base; - } - - void shrinkStep(Node& top, Node& middle, Node& bottom, - typename Graph::template NodeMap& ear, - UFE& blossom, NV& rep, EFE& tree, std::queue& Q) { - //We traverse a blossom and update everything. - - ear.set(top,bottom); - Node t=top; - while ( t!=middle ) { - Node u=_mate[t]; - t=ear[u]; - ear.set(t,u); - } - bottom=_mate[middle]; - position.set(bottom,D); - Q.push(bottom); - top=ear[bottom]; - Node oldmiddle=middle; - middle=rep[blossom.find(top)]; - tree.erase(bottom); - tree.erase(oldmiddle); - blossom.insert(bottom); - blossom.join(bottom, oldmiddle); - blossom.join(top, oldmiddle); - } - - - - bool growOrAugment(Node& y, Node& x, typename Graph::template - NodeMap& ear, UFE& blossom, NV& rep, EFE& tree, - std::queue& Q) { - //x is in a blossom in the tree, y is outside. If y is covered by - //the matching we grow, otherwise we augment. In this case we - //return 1. - - if ( _mate[y]!=INVALID ) { //grow - ear.set(y,x); - Node w=_mate[y]; - rep[blossom.insert(w)] = w; - position.set(y,A); - position.set(w,D); - int t = tree.find(rep[blossom.find(x)]); - tree.insert(y,t); - tree.insert(w,t); - Q.push(w); - } else { //augment - augment(x, ear, blossom, rep, tree); - _mate.set(x,y); - _mate.set(y,x); - return true; - } - return false; - } - - void augment(Node x, typename Graph::template NodeMap& ear, - UFE& blossom, NV& rep, EFE& tree) { - Node v=_mate[x]; - while ( v!=INVALID ) { - - Node u=ear[v]; - _mate.set(v,u); - Node tmp=v; - v=_mate[u]; - _mate.set(u,tmp); - } - int y = tree.find(rep[blossom.find(x)]); - for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) { - if ( position[tit] == D ) { - int b = blossom.find(tit); - for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) { - position.set(bit, C); - } - blossom.eraseClass(b); - } else position.set(tit, C); - } - tree.eraseClass(y); - - } + /// @} }; @@ -627,25 +617,28 @@ /// \f$O(nm\log(n))\f$ time complexity. /// /// The maximum weighted matching problem is to find undirected - /// arcs in the digraph with maximum overall weight and no two of - /// them shares their endpoints. The problem can be formulated with - /// the next linear program: + /// edges in the graph with maximum overall weight and no two of + /// them shares their ends. The problem can be formulated with the + /// following linear program. /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f] - ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f] + /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} + \quad \forall B\in\mathcal{O}\f] */ /// \f[x_e \ge 0\quad \forall e\in E\f] /// \f[\max \sum_{e\in E}x_ew_e\f] - /// where \f$\delta(X)\f$ is the set of arcs incident to a node in - /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in - /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of - /// the nodes. + /// where \f$\delta(X)\f$ is the set of edges incident to a node in + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality + /// subsets of the nodes. /// /// The algorithm calculates an optimal matching and a proof of the /// optimality. The solution of the dual problem can be used to check - /// the result of the algorithm. The dual linear problem is the next: - /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f] + /// the result of the algorithm. The dual linear problem is the + /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)} + z_B \ge w_{uv} \quad \forall uv\in E\f] */ /// \f[y_u \ge 0 \quad \forall u \in V\f] /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] - /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f] + /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} + \frac{\vert B \vert - 1}{2}z_B\f] */ /// /// The algorithm can be executed with \c run() or the \c init() and /// then the \c start() member functions. After it the matching can @@ -705,16 +698,13 @@ int _node_num; int _blossom_num; - typedef typename Graph::template NodeMap NodeIntMap; - typedef typename Graph::template ArcMap ArcIntMap; - typedef typename Graph::template EdgeMap EdgeIntMap; typedef RangeMap IntIntMap; enum Status { EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2 }; - typedef HeapUnionFind BlossomSet; + typedef HeapUnionFind BlossomSet; struct BlossomData { int tree; Status status; @@ -723,21 +713,21 @@ Node base; }; - NodeIntMap *_blossom_index; + IntNodeMap *_blossom_index; BlossomSet *_blossom_set; RangeMap* _blossom_data; - NodeIntMap *_node_index; - ArcIntMap *_node_heap_index; + IntNodeMap *_node_index; + IntArcMap *_node_heap_index; struct NodeData { - NodeData(ArcIntMap& node_heap_index) + NodeData(IntArcMap& node_heap_index) : heap(node_heap_index) {} int blossom; Value pot; - BinHeap heap; + BinHeap heap; std::map heap_index; int tree; @@ -750,14 +740,14 @@ IntIntMap *_tree_set_index; TreeSet *_tree_set; - NodeIntMap *_delta1_index; - BinHeap *_delta1; + IntNodeMap *_delta1_index; + BinHeap *_delta1; IntIntMap *_delta2_index; BinHeap *_delta2; - EdgeIntMap *_delta3_index; - BinHeap *_delta3; + IntEdgeMap *_delta3_index; + BinHeap *_delta3; IntIntMap *_delta4_index; BinHeap *_delta4; @@ -775,14 +765,14 @@ _node_potential = new NodePotential(_graph); } if (!_blossom_set) { - _blossom_index = new NodeIntMap(_graph); + _blossom_index = new IntNodeMap(_graph); _blossom_set = new BlossomSet(*_blossom_index); _blossom_data = new RangeMap(_blossom_num); } if (!_node_index) { - _node_index = new NodeIntMap(_graph); - _node_heap_index = new ArcIntMap(_graph); + _node_index = new IntNodeMap(_graph); + _node_heap_index = new IntArcMap(_graph); _node_data = new RangeMap(_node_num, NodeData(*_node_heap_index)); } @@ -792,16 +782,16 @@ _tree_set = new TreeSet(*_tree_set_index); } if (!_delta1) { - _delta1_index = new NodeIntMap(_graph); - _delta1 = new BinHeap(*_delta1_index); + _delta1_index = new IntNodeMap(_graph); + _delta1 = new BinHeap(*_delta1_index); } if (!_delta2) { _delta2_index = new IntIntMap(_blossom_num); _delta2 = new BinHeap(*_delta2_index); } if (!_delta3) { - _delta3_index = new EdgeIntMap(_graph); - _delta3 = new BinHeap(*_delta3_index); + _delta3_index = new IntEdgeMap(_graph); + _delta3 = new BinHeap(*_delta3_index); } if (!_delta4) { _delta4_index = new IntIntMap(_blossom_num); @@ -1266,10 +1256,10 @@ } - void augmentOnArc(const Edge& arc) { - - int left = _blossom_set->find(_graph.u(arc)); - int right = _blossom_set->find(_graph.v(arc)); + void augmentOnEdge(const Edge& edge) { + + int left = _blossom_set->find(_graph.u(edge)); + int right = _blossom_set->find(_graph.v(edge)); if ((*_blossom_data)[left].status == EVEN) { int left_tree = _tree_set->find(left); @@ -1289,8 +1279,8 @@ unmatchedToMatched(right); } - (*_blossom_data)[left].next = _graph.direct(arc, true); - (*_blossom_data)[right].next = _graph.direct(arc, false); + (*_blossom_data)[left].next = _graph.direct(edge, true); + (*_blossom_data)[right].next = _graph.direct(edge, false); } void extendOnArc(const Arc& arc) { @@ -1310,7 +1300,7 @@ matchedToEven(even, tree); } - void shrinkOnArc(const Edge& edge, int tree) { + void shrinkOnEdge(const Edge& edge, int tree) { int nca = -1; std::vector left_path, right_path; @@ -1652,7 +1642,7 @@ createStructures(); for (ArcIt e(_graph); e != INVALID; ++e) { - _node_heap_index->set(e, BinHeap::PRE_HEAP); + _node_heap_index->set(e, BinHeap::PRE_HEAP); } for (NodeIt n(_graph); n != INVALID; ++n) { _delta1_index->set(n, _delta1->PRE_HEAP); @@ -1769,9 +1759,9 @@ } if (left_tree == right_tree) { - shrinkOnArc(e, left_tree); + shrinkOnEdge(e, left_tree); } else { - augmentOnArc(e); + augmentOnEdge(e); unmatched -= 2; } } @@ -1818,11 +1808,24 @@ return sum /= 2; } - /// \brief Returns true when the arc is in the matching. + /// \brief Returns the cardinality of the matching. /// - /// Returns true when the arc is in the matching. - bool matching(const Edge& arc) const { - return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true); + /// Returns the cardinality of the matching. + int matchingSize() const { + int num = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_matching)[n] != INVALID) { + ++num; + } + } + return num /= 2; + } + + /// \brief Returns true when the edge is in the matching. + /// + /// Returns true when the edge is in the matching. + bool matching(const Edge& edge) const { + return edge == (*_matching)[_graph.u(edge)]; } /// \brief Returns the incident matching arc. @@ -1913,16 +1916,11 @@ _last = _algorithm->_blossom_potential[variable].end; } - /// \brief Invalid constructor. - /// - /// Invalid constructor. - BlossomIt(Invalid) : _index(-1) {} - /// \brief Conversion to node. /// /// Conversion to node. operator Node() const { - return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID; + return _algorithm->_blossom_node_list[_index]; } /// \brief Increment operator. @@ -1930,18 +1928,18 @@ /// Increment operator. BlossomIt& operator++() { ++_index; - if (_index == _last) { - _index = -1; - } return *this; } - bool operator==(const BlossomIt& it) const { - return _index == it._index; - } - bool operator!=(const BlossomIt& it) const { - return _index != it._index; - } + /// \brief Validity checking + /// + /// Checks whether the iterator is invalid. + bool operator==(Invalid) const { return _index == _last; } + + /// \brief Validity checking + /// + /// Checks whether the iterator is valid. + bool operator!=(Invalid) const { return _index != _last; } private: const MaxWeightedMatching* _algorithm; @@ -1958,29 +1956,32 @@ /// \brief Weighted perfect matching in general graphs /// /// This class provides an efficient implementation of Edmond's - /// maximum weighted perfecr matching algorithm. The implementation + /// maximum weighted perfect matching algorithm. The implementation /// is based on extensive use of priority queues and provides /// \f$O(nm\log(n))\f$ time complexity. /// /// The maximum weighted matching problem is to find undirected - /// arcs in the digraph with maximum overall weight and no two of - /// them shares their endpoints and covers all nodes. The problem - /// can be formulated with the next linear program: + /// edges in the graph with maximum overall weight and no two of + /// them shares their ends and covers all nodes. The problem can be + /// formulated with the following linear program. /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f] - ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f] + /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} + \quad \forall B\in\mathcal{O}\f] */ /// \f[x_e \ge 0\quad \forall e\in E\f] /// \f[\max \sum_{e\in E}x_ew_e\f] - /// where \f$\delta(X)\f$ is the set of arcs incident to a node in - /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in - /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of - /// the nodes. + /// where \f$\delta(X)\f$ is the set of edges incident to a node in + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality + /// subsets of the nodes. /// /// The algorithm calculates an optimal matching and a proof of the /// optimality. The solution of the dual problem can be used to check - /// the result of the algorithm. The dual linear problem is the next: - /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f] + /// the result of the algorithm. The dual linear problem is the + /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge + w_{uv} \quad \forall uv\in E\f] */ /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] - /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f] + /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} + \frac{\vert B \vert - 1}{2}z_B\f] */ /// /// The algorithm can be executed with \c run() or the \c init() and /// then the \c start() member functions. After it the matching can @@ -2040,16 +2041,13 @@ int _node_num; int _blossom_num; - typedef typename Graph::template NodeMap NodeIntMap; - typedef typename Graph::template ArcMap ArcIntMap; - typedef typename Graph::template EdgeMap EdgeIntMap; typedef RangeMap IntIntMap; enum Status { EVEN = -1, MATCHED = 0, ODD = 1 }; - typedef HeapUnionFind BlossomSet; + typedef HeapUnionFind BlossomSet; struct BlossomData { int tree; Status status; @@ -2057,21 +2055,21 @@ Value pot, offset; }; - NodeIntMap *_blossom_index; + IntNodeMap *_blossom_index; BlossomSet *_blossom_set; RangeMap* _blossom_data; - NodeIntMap *_node_index; - ArcIntMap *_node_heap_index; + IntNodeMap *_node_index; + IntArcMap *_node_heap_index; struct NodeData { - NodeData(ArcIntMap& node_heap_index) + NodeData(IntArcMap& node_heap_index) : heap(node_heap_index) {} int blossom; Value pot; - BinHeap heap; + BinHeap heap; std::map heap_index; int tree; @@ -2087,8 +2085,8 @@ IntIntMap *_delta2_index; BinHeap *_delta2; - EdgeIntMap *_delta3_index; - BinHeap *_delta3; + IntEdgeMap *_delta3_index; + BinHeap *_delta3; IntIntMap *_delta4_index; BinHeap *_delta4; @@ -2106,16 +2104,16 @@ _node_potential = new NodePotential(_graph); } if (!_blossom_set) { - _blossom_index = new NodeIntMap(_graph); + _blossom_index = new IntNodeMap(_graph); _blossom_set = new BlossomSet(*_blossom_index); _blossom_data = new RangeMap(_blossom_num); } if (!_node_index) { - _node_index = new NodeIntMap(_graph); - _node_heap_index = new ArcIntMap(_graph); + _node_index = new IntNodeMap(_graph); + _node_heap_index = new IntArcMap(_graph); _node_data = new RangeMap(_node_num, - NodeData(*_node_heap_index)); + NodeData(*_node_heap_index)); } if (!_tree_set) { @@ -2127,8 +2125,8 @@ _delta2 = new BinHeap(*_delta2_index); } if (!_delta3) { - _delta3_index = new EdgeIntMap(_graph); - _delta3 = new BinHeap(*_delta3_index); + _delta3_index = new IntEdgeMap(_graph); + _delta3 = new BinHeap(*_delta3_index); } if (!_delta4) { _delta4_index = new IntIntMap(_blossom_num); @@ -2461,10 +2459,10 @@ _tree_set->eraseClass(tree); } - void augmentOnArc(const Edge& arc) { - - int left = _blossom_set->find(_graph.u(arc)); - int right = _blossom_set->find(_graph.v(arc)); + void augmentOnEdge(const Edge& edge) { + + int left = _blossom_set->find(_graph.u(edge)); + int right = _blossom_set->find(_graph.v(edge)); int left_tree = _tree_set->find(left); alternatePath(left, left_tree); @@ -2474,8 +2472,8 @@ alternatePath(right, right_tree); destroyTree(right_tree); - (*_blossom_data)[left].next = _graph.direct(arc, true); - (*_blossom_data)[right].next = _graph.direct(arc, false); + (*_blossom_data)[left].next = _graph.direct(edge, true); + (*_blossom_data)[right].next = _graph.direct(edge, false); } void extendOnArc(const Arc& arc) { @@ -2495,7 +2493,7 @@ matchedToEven(even, tree); } - void shrinkOnArc(const Edge& edge, int tree) { + void shrinkOnEdge(const Edge& edge, int tree) { int nca = -1; std::vector left_path, right_path; @@ -2831,7 +2829,7 @@ createStructures(); for (ArcIt e(_graph); e != INVALID; ++e) { - _node_heap_index->set(e, BinHeap::PRE_HEAP); + _node_heap_index->set(e, BinHeap::PRE_HEAP); } for (EdgeIt e(_graph); e != INVALID; ++e) { _delta3_index->set(e, _delta3->PRE_HEAP); @@ -2924,9 +2922,9 @@ int right_tree = _tree_set->find(right_blossom); if (left_tree == right_tree) { - shrinkOnArc(e, left_tree); + shrinkOnEdge(e, left_tree); } else { - augmentOnArc(e); + augmentOnEdge(e); unmatched -= 2; } } @@ -2974,16 +2972,16 @@ return sum /= 2; } - /// \brief Returns true when the arc is in the matching. + /// \brief Returns true when the edge is in the matching. /// - /// Returns true when the arc is in the matching. - bool matching(const Edge& arc) const { - return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true); + /// Returns true when the edge is in the matching. + bool matching(const Edge& edge) const { + return static_cast((*_matching)[_graph.u(edge)]) == edge; } - /// \brief Returns the incident matching arc. + /// \brief Returns the incident matching edge. /// - /// Returns the incident matching arc from given node. + /// Returns the incident matching arc from given edge. Arc matching(const Node& node) const { return (*_matching)[node]; } @@ -3066,16 +3064,11 @@ _last = _algorithm->_blossom_potential[variable].end; } - /// \brief Invalid constructor. - /// - /// Invalid constructor. - BlossomIt(Invalid) : _index(-1) {} - /// \brief Conversion to node. /// /// Conversion to node. operator Node() const { - return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID; + return _algorithm->_blossom_node_list[_index]; } /// \brief Increment operator. @@ -3083,18 +3076,18 @@ /// Increment operator. BlossomIt& operator++() { ++_index; - if (_index == _last) { - _index = -1; - } return *this; } - bool operator==(const BlossomIt& it) const { - return _index == it._index; - } - bool operator!=(const BlossomIt& it) const { - return _index != it._index; - } + /// \brief Validity checking + /// + /// Checks whether the iterator is invalid. + bool operator==(Invalid) const { return _index == _last; } + + /// \brief Validity checking + /// + /// Checks whether the iterator is valid. + bool operator!=(Invalid) const { return _index != _last; } private: const MaxWeightedPerfectMatching* _algorithm; diff -r 64ad48007fb2 -r 91d63b8b1a4c test/CMakeLists.txt --- a/test/CMakeLists.txt Mon Oct 13 13:56:00 2008 +0200 +++ b/test/CMakeLists.txt Mon Oct 13 14:00:11 2008 +0200 @@ -17,7 +17,6 @@ kruskal_test maps_test max_matching_test - max_weighted_matching_test random_test path_test time_measure_test diff -r 64ad48007fb2 -r 91d63b8b1a4c test/Makefile.am --- a/test/Makefile.am Mon Oct 13 13:56:00 2008 +0200 +++ b/test/Makefile.am Mon Oct 13 14:00:11 2008 +0200 @@ -20,7 +20,6 @@ test/kruskal_test \ test/maps_test \ test/max_matching_test \ - test/max_weighted_matching_test \ test/random_test \ test/path_test \ test/test_tools_fail \ @@ -45,7 +44,6 @@ test_kruskal_test_SOURCES = test/kruskal_test.cc test_maps_test_SOURCES = test/maps_test.cc test_max_matching_test_SOURCES = test/max_matching_test.cc -test_max_weighted_matching_test_SOURCES = test/max_weighted_matching_test.cc test_path_test_SOURCES = test/path_test.cc test_random_test_SOURCES = test/random_test.cc test_test_tools_fail_SOURCES = test/test_tools_fail.cc diff -r 64ad48007fb2 -r 91d63b8b1a4c test/max_matching_test.cc --- a/test/max_matching_test.cc Mon Oct 13 13:56:00 2008 +0200 +++ b/test/max_matching_test.cc Mon Oct 13 14:00:11 2008 +0200 @@ -17,165 +17,294 @@ */ #include +#include #include #include #include #include +#include +#include +#include + #include "test_tools.h" -#include -#include using namespace std; using namespace lemon; +GRAPH_TYPEDEFS(SmartGraph); + + +const int lgfn = 3; +const std::string lgf[lgfn] = { + "@nodes\n" + "label\n" + "0\n" + "1\n" + "2\n" + "3\n" + "4\n" + "5\n" + "6\n" + "7\n" + "@edges\n" + " label weight\n" + "7 4 0 984\n" + "0 7 1 73\n" + "7 1 2 204\n" + "2 3 3 583\n" + "2 7 4 565\n" + "2 1 5 582\n" + "0 4 6 551\n" + "2 5 7 385\n" + "1 5 8 561\n" + "5 3 9 484\n" + "7 5 10 904\n" + "3 6 11 47\n" + "7 6 12 888\n" + "3 0 13 747\n" + "6 1 14 310\n", + + "@nodes\n" + "label\n" + "0\n" + "1\n" + "2\n" + "3\n" + "4\n" + "5\n" + "6\n" + "7\n" + "@edges\n" + " label weight\n" + "2 5 0 710\n" + "0 5 1 241\n" + "2 4 2 856\n" + "2 6 3 762\n" + "4 1 4 747\n" + "6 1 5 962\n" + "4 7 6 723\n" + "1 7 7 661\n" + "2 3 8 376\n" + "1 0 9 416\n" + "6 7 10 391\n", + + "@nodes\n" + "label\n" + "0\n" + "1\n" + "2\n" + "3\n" + "4\n" + "5\n" + "6\n" + "7\n" + "@edges\n" + " label weight\n" + "6 2 0 553\n" + "0 7 1 653\n" + "6 3 2 22\n" + "4 7 3 846\n" + "7 2 4 981\n" + "7 6 5 250\n" + "5 2 6 539\n", +}; + +void checkMatching(const SmartGraph& graph, + const MaxMatching& mm) { + int num = 0; + + IntNodeMap comp_index(graph); + UnionFind comp(comp_index); + + int barrier_num = 0; + + for (NodeIt n(graph); n != INVALID; ++n) { + check(mm.decomposition(n) == MaxMatching::EVEN || + mm.matching(n) != INVALID, "Wrong Gallai-Edmonds decomposition"); + if (mm.decomposition(n) == MaxMatching::ODD) { + ++barrier_num; + } else { + comp.insert(n); + } + } + + for (EdgeIt e(graph); e != INVALID; ++e) { + if (mm.matching(e)) { + check(e == mm.matching(graph.u(e)), "Wrong matching"); + check(e == mm.matching(graph.v(e)), "Wrong matching"); + ++num; + } + check(mm.decomposition(graph.u(e)) != MaxMatching::EVEN || + mm.decomposition(graph.v(e)) != MaxMatching::MATCHED, + "Wrong Gallai-Edmonds decomposition"); + + check(mm.decomposition(graph.v(e)) != MaxMatching::EVEN || + mm.decomposition(graph.u(e)) != MaxMatching::MATCHED, + "Wrong Gallai-Edmonds decomposition"); + + if (mm.decomposition(graph.u(e)) != MaxMatching::ODD && + mm.decomposition(graph.v(e)) != MaxMatching::ODD) { + comp.join(graph.u(e), graph.v(e)); + } + } + + std::set comp_root; + int odd_comp_num = 0; + for (NodeIt n(graph); n != INVALID; ++n) { + if (mm.decomposition(n) != MaxMatching::ODD) { + int root = comp.find(n); + if (comp_root.find(root) == comp_root.end()) { + comp_root.insert(root); + if (comp.size(n) % 2 == 1) { + ++odd_comp_num; + } + } + } + } + + check(mm.matchingSize() == num, "Wrong matching"); + check(2 * num == countNodes(graph) - (odd_comp_num - barrier_num), + "Wrong matching"); + return; +} + +void checkWeightedMatching(const SmartGraph& graph, + const SmartGraph::EdgeMap& weight, + const MaxWeightedMatching& mwm) { + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) { + if (graph.u(e) == graph.v(e)) continue; + int rw = mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e)); + + for (int i = 0; i < mwm.blossomNum(); ++i) { + bool s = false, t = false; + for (MaxWeightedMatching::BlossomIt n(mwm, i); + n != INVALID; ++n) { + if (graph.u(e) == n) s = true; + if (graph.v(e) == n) t = true; + } + if (s == true && t == true) { + rw += mwm.blossomValue(i); + } + } + rw -= weight[e] * mwm.dualScale; + + check(rw >= 0, "Negative reduced weight"); + check(rw == 0 || !mwm.matching(e), + "Non-zero reduced weight on matching edge"); + } + + int pv = 0; + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + if (mwm.matching(n) != INVALID) { + check(mwm.nodeValue(n) >= 0, "Invalid node value"); + pv += weight[mwm.matching(n)]; + SmartGraph::Node o = graph.target(mwm.matching(n)); + check(mwm.mate(n) == o, "Invalid matching"); + check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)), + "Invalid matching"); + } else { + check(mwm.mate(n) == INVALID, "Invalid matching"); + check(mwm.nodeValue(n) == 0, "Invalid matching"); + } + } + + int dv = 0; + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + dv += mwm.nodeValue(n); + } + + for (int i = 0; i < mwm.blossomNum(); ++i) { + check(mwm.blossomValue(i) >= 0, "Invalid blossom value"); + check(mwm.blossomSize(i) % 2 == 1, "Even blossom size"); + dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2); + } + + check(pv * mwm.dualScale == dv * 2, "Wrong duality"); + + return; +} + +void checkWeightedPerfectMatching(const SmartGraph& graph, + const SmartGraph::EdgeMap& weight, + const MaxWeightedPerfectMatching& mwpm) { + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) { + if (graph.u(e) == graph.v(e)) continue; + int rw = mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e)); + + for (int i = 0; i < mwpm.blossomNum(); ++i) { + bool s = false, t = false; + for (MaxWeightedPerfectMatching::BlossomIt n(mwpm, i); + n != INVALID; ++n) { + if (graph.u(e) == n) s = true; + if (graph.v(e) == n) t = true; + } + if (s == true && t == true) { + rw += mwpm.blossomValue(i); + } + } + rw -= weight[e] * mwpm.dualScale; + + check(rw >= 0, "Negative reduced weight"); + check(rw == 0 || !mwpm.matching(e), + "Non-zero reduced weight on matching edge"); + } + + int pv = 0; + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + check(mwpm.matching(n) != INVALID, "Non perfect"); + pv += weight[mwpm.matching(n)]; + SmartGraph::Node o = graph.target(mwpm.matching(n)); + check(mwpm.mate(n) == o, "Invalid matching"); + check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)), + "Invalid matching"); + } + + int dv = 0; + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + dv += mwpm.nodeValue(n); + } + + for (int i = 0; i < mwpm.blossomNum(); ++i) { + check(mwpm.blossomValue(i) >= 0, "Invalid blossom value"); + check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size"); + dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2); + } + + check(pv * mwpm.dualScale == dv * 2, "Wrong duality"); + + return; +} + + int main() { - typedef ListGraph Graph; + for (int i = 0; i < lgfn; ++i) { + SmartGraph graph; + SmartGraph::EdgeMap weight(graph); - GRAPH_TYPEDEFS(Graph); + istringstream lgfs(lgf[i]); + graphReader(graph, lgfs). + edgeMap("weight", weight).run(); - Graph g; - g.clear(); + MaxMatching mm(graph); + mm.run(); + checkMatching(graph, mm); - std::vector nodes; - for (int i=0; i<13; ++i) - nodes.push_back(g.addNode()); + MaxWeightedMatching mwm(graph, weight); + mwm.run(); + checkWeightedMatching(graph, weight, mwm); - g.addEdge(nodes[0], nodes[0]); - g.addEdge(nodes[6], nodes[10]); - g.addEdge(nodes[5], nodes[10]); - g.addEdge(nodes[4], nodes[10]); - g.addEdge(nodes[3], nodes[11]); - g.addEdge(nodes[1], nodes[6]); - g.addEdge(nodes[4], nodes[7]); - g.addEdge(nodes[1], nodes[8]); - g.addEdge(nodes[0], nodes[8]); - g.addEdge(nodes[3], nodes[12]); - g.addEdge(nodes[6], nodes[9]); - g.addEdge(nodes[9], nodes[11]); - g.addEdge(nodes[2], nodes[10]); - g.addEdge(nodes[10], nodes[8]); - g.addEdge(nodes[5], nodes[8]); - g.addEdge(nodes[6], nodes[3]); - g.addEdge(nodes[0], nodes[5]); - g.addEdge(nodes[6], nodes[12]); + MaxWeightedPerfectMatching mwpm(graph, weight); + bool perfect = mwpm.run(); - MaxMatching max_matching(g); - max_matching.init(); - max_matching.startDense(); + check(perfect == (mm.matchingSize() * 2 == countNodes(graph)), + "Perfect matching found"); - int s=0; - Graph::NodeMap mate(g,INVALID); - max_matching.mateMap(mate); - for(NodeIt v(g); v!=INVALID; ++v) { - if ( mate[v]!=INVALID ) ++s; - } - int size=int(s/2); //size will be used as the size of a maxmatching - - for(NodeIt v(g); v!=INVALID; ++v) { - max_matching.mate(v); - } - - check ( size == max_matching.size(), "mate() returns a different size matching than max_matching.size()" ); - - Graph::NodeMap::DecompType> pos0(g); - max_matching.decomposition(pos0); - - max_matching.init(); - max_matching.startSparse(); - s=0; - max_matching.mateMap(mate); - for(NodeIt v(g); v!=INVALID; ++v) { - if ( mate[v]!=INVALID ) ++s; - } - check ( int(s/2) == size, "The size does not equal!" ); - - Graph::NodeMap::DecompType> pos1(g); - max_matching.decomposition(pos1); - - max_matching.run(); - s=0; - max_matching.mateMap(mate); - for(NodeIt v(g); v!=INVALID; ++v) { - if ( mate[v]!=INVALID ) ++s; - } - check ( int(s/2) == size, "The size does not equal!" ); - - Graph::NodeMap::DecompType> pos2(g); - max_matching.decomposition(pos2); - - max_matching.run(); - s=0; - max_matching.mateMap(mate); - for(NodeIt v(g); v!=INVALID; ++v) { - if ( mate[v]!=INVALID ) ++s; - } - check ( int(s/2) == size, "The size does not equal!" ); - - Graph::NodeMap::DecompType> pos(g); - max_matching.decomposition(pos); - - bool ismatching=true; - for(NodeIt v(g); v!=INVALID; ++v) { - if ( mate[v]!=INVALID ) { - Node u=mate[v]; - if (mate[u]!=v) ismatching=false; + if (perfect) { + checkWeightedPerfectMatching(graph, weight, mwpm); } } - check ( ismatching, "It is not a matching!" ); - - bool coincide=true; - for(NodeIt v(g); v!=INVALID; ++v) { - if ( pos0[v] != pos1[v] || pos1[v]!=pos2[v] || pos2[v]!=pos[v] ) { - coincide=false; - } - } - check ( coincide, "The decompositions do not coincide! " ); - - bool noarc=true; - for(EdgeIt e(g); e!=INVALID; ++e) { - if ( (pos[g.v(e)]==max_matching.C && - pos[g.u(e)]==max_matching.D) || - (pos[g.v(e)]==max_matching.D && - pos[g.u(e)]==max_matching.C) ) - noarc=false; - } - check ( noarc, "There are arcs between D and C!" ); - - bool oddcomp=true; - Graph::NodeMap todo(g,true); - int num_comp=0; - for(NodeIt v(g); v!=INVALID; ++v) { - if ( pos[v]==max_matching.D && todo[v] ) { - int comp_size=1; - ++num_comp; - std::queue Q; - Q.push(v); - todo.set(v,false); - while (!Q.empty()) { - Node w=Q.front(); - Q.pop(); - for(IncEdgeIt e(g,w); e!=INVALID; ++e) { - Node u=g.runningNode(e); - if ( pos[u]==max_matching.D && todo[u] ) { - ++comp_size; - Q.push(u); - todo.set(u,false); - } - } - } - if ( !(comp_size % 2) ) oddcomp=false; - } - } - check ( oddcomp, "A component of g[D] is not odd." ); - - int barrier=0; - for(NodeIt v(g); v!=INVALID; ++v) { - if ( pos[v]==max_matching.A ) ++barrier; - } - int expected_size=int( countNodes(g)-num_comp+barrier)/2; - check ( size==expected_size, "The size of the matching is wrong." ); return 0; } diff -r 64ad48007fb2 -r 91d63b8b1a4c test/max_weighted_matching_test.cc --- a/test/max_weighted_matching_test.cc Mon Oct 13 13:56:00 2008 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,250 +0,0 @@ -/* -*- mode: C++; indent-tabs-mode: nil; -*- - * - * This file is a part of LEMON, a generic C++ optimization library. - * - * Copyright (C) 2003-2008 - * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport - * (Egervary Research Group on Combinatorial Optimization, EGRES). - * - * Permission to use, modify and distribute this software is granted - * provided that this copyright notice appears in all copies. For - * precise terms see the accompanying LICENSE file. - * - * This software is provided "AS IS" with no warranty of any kind, - * express or implied, and with no claim as to its suitability for any - * purpose. - * - */ - -#include -#include -#include -#include -#include -#include - -#include "test_tools.h" -#include -#include -#include - -using namespace std; -using namespace lemon; - -GRAPH_TYPEDEFS(SmartGraph); - -const int lgfn = 3; -const std::string lgf[lgfn] = { - "@nodes\n" - "label\n" - "0\n" - "1\n" - "2\n" - "3\n" - "4\n" - "5\n" - "6\n" - "7\n" - "@edges\n" - "label weight\n" - "7 4 0 984\n" - "0 7 1 73\n" - "7 1 2 204\n" - "2 3 3 583\n" - "2 7 4 565\n" - "2 1 5 582\n" - "0 4 6 551\n" - "2 5 7 385\n" - "1 5 8 561\n" - "5 3 9 484\n" - "7 5 10 904\n" - "3 6 11 47\n" - "7 6 12 888\n" - "3 0 13 747\n" - "6 1 14 310\n", - - "@nodes\n" - "label\n" - "0\n" - "1\n" - "2\n" - "3\n" - "4\n" - "5\n" - "6\n" - "7\n" - "@edges\n" - "label weight\n" - "2 5 0 710\n" - "0 5 1 241\n" - "2 4 2 856\n" - "2 6 3 762\n" - "4 1 4 747\n" - "6 1 5 962\n" - "4 7 6 723\n" - "1 7 7 661\n" - "2 3 8 376\n" - "1 0 9 416\n" - "6 7 10 391\n", - - "@nodes\n" - "label\n" - "0\n" - "1\n" - "2\n" - "3\n" - "4\n" - "5\n" - "6\n" - "7\n" - "@edges\n" - "label weight\n" - "6 2 0 553\n" - "0 7 1 653\n" - "6 3 2 22\n" - "4 7 3 846\n" - "7 2 4 981\n" - "7 6 5 250\n" - "5 2 6 539\n" -}; - -void checkMatching(const SmartGraph& graph, - const SmartGraph::EdgeMap& weight, - const MaxWeightedMatching& mwm) { - for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) { - if (graph.u(e) == graph.v(e)) continue; - int rw = - mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e)); - - for (int i = 0; i < mwm.blossomNum(); ++i) { - bool s = false, t = false; - for (MaxWeightedMatching::BlossomIt n(mwm, i); - n != INVALID; ++n) { - if (graph.u(e) == n) s = true; - if (graph.v(e) == n) t = true; - } - if (s == true && t == true) { - rw += mwm.blossomValue(i); - } - } - rw -= weight[e] * mwm.dualScale; - - check(rw >= 0, "Negative reduced weight"); - check(rw == 0 || !mwm.matching(e), - "Non-zero reduced weight on matching arc"); - } - - int pv = 0; - for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { - if (mwm.matching(n) != INVALID) { - check(mwm.nodeValue(n) >= 0, "Invalid node value"); - pv += weight[mwm.matching(n)]; - SmartGraph::Node o = graph.target(mwm.matching(n)); - check(mwm.mate(n) == o, "Invalid matching"); - check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)), - "Invalid matching"); - } else { - check(mwm.mate(n) == INVALID, "Invalid matching"); - check(mwm.nodeValue(n) == 0, "Invalid matching"); - } - } - - int dv = 0; - for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { - dv += mwm.nodeValue(n); - } - - for (int i = 0; i < mwm.blossomNum(); ++i) { - check(mwm.blossomValue(i) >= 0, "Invalid blossom value"); - check(mwm.blossomSize(i) % 2 == 1, "Even blossom size"); - dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2); - } - - check(pv * mwm.dualScale == dv * 2, "Wrong duality"); - - return; -} - -void checkPerfectMatching(const SmartGraph& graph, - const SmartGraph::EdgeMap& weight, - const MaxWeightedPerfectMatching& mwpm) { - for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) { - if (graph.u(e) == graph.v(e)) continue; - int rw = - mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e)); - - for (int i = 0; i < mwpm.blossomNum(); ++i) { - bool s = false, t = false; - for (MaxWeightedPerfectMatching::BlossomIt n(mwpm, i); - n != INVALID; ++n) { - if (graph.u(e) == n) s = true; - if (graph.v(e) == n) t = true; - } - if (s == true && t == true) { - rw += mwpm.blossomValue(i); - } - } - rw -= weight[e] * mwpm.dualScale; - - check(rw >= 0, "Negative reduced weight"); - check(rw == 0 || !mwpm.matching(e), - "Non-zero reduced weight on matching arc"); - } - - int pv = 0; - for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { - check(mwpm.matching(n) != INVALID, "Non perfect"); - pv += weight[mwpm.matching(n)]; - SmartGraph::Node o = graph.target(mwpm.matching(n)); - check(mwpm.mate(n) == o, "Invalid matching"); - check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)), - "Invalid matching"); - } - - int dv = 0; - for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { - dv += mwpm.nodeValue(n); - } - - for (int i = 0; i < mwpm.blossomNum(); ++i) { - check(mwpm.blossomValue(i) >= 0, "Invalid blossom value"); - check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size"); - dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2); - } - - check(pv * mwpm.dualScale == dv * 2, "Wrong duality"); - - return; -} - - -int main() { - - for (int i = 0; i < lgfn; ++i) { - SmartGraph graph; - SmartGraph::EdgeMap weight(graph); - - istringstream lgfs(lgf[i]); - graphReader(graph, lgfs). - edgeMap("weight", weight).run(); - - MaxWeightedMatching mwm(graph, weight); - mwm.run(); - checkMatching(graph, weight, mwm); - - MaxMatching mm(graph); - mm.run(); - - MaxWeightedPerfectMatching mwpm(graph, weight); - bool perfect = mwpm.run(); - - check(perfect == (mm.size() * 2 == countNodes(graph)), - "Perfect matching found"); - - if (perfect) { - checkPerfectMatching(graph, weight, mwpm); - } - } - - return 0; -}