diff --git a/doc/groups.dox b/doc/groups.dox --- a/doc/groups.dox +++ b/doc/groups.dox @@ -435,9 +435,10 @@ @ingroup algs \brief Algorithms for finding matchings in graphs and bipartite graphs. -This group contains algorithm objects and functions to calculate +This group contains the algorithms for calculating matchings in graphs and bipartite graphs. The general matching problem is -finding a subset of the arcs which does not shares common endpoints. +finding a subset of the edges for which each node has at most one incident +edge. There are several different algorithms for calculate matchings in graphs. The matching problems in bipartite graphs are generally diff --git a/lemon/Makefile.am b/lemon/Makefile.am --- a/lemon/Makefile.am +++ b/lemon/Makefile.am @@ -89,8 +89,8 @@ lemon/lp_skeleton.h \ lemon/list_graph.h \ lemon/maps.h \ + lemon/matching.h \ lemon/math.h \ - lemon/max_matching.h \ lemon/min_cost_arborescence.h \ lemon/nauty_reader.h \ lemon/path.h \ diff --git a/lemon/euler.h b/lemon/euler.h --- a/lemon/euler.h +++ b/lemon/euler.h @@ -26,33 +26,31 @@ /// \ingroup graph_properties /// \file -/// \brief Euler tour +/// \brief Euler tour iterators and a function for checking the \e Eulerian +/// property. /// -///This file provides an Euler tour iterator and ways to check -///if a digraph is euler. - +///This file provides Euler tour iterators and a function to check +///if a (di)graph is \e Eulerian. namespace lemon { - ///Euler iterator for digraphs. + ///Euler tour iterator for digraphs. - /// \ingroup graph_properties - ///This iterator converts to the \c Arc type of the digraph and using - ///operator ++, it provides an Euler tour of a \e directed - ///graph (if there exists). + /// \ingroup graph_prop + ///This iterator provides an Euler tour (Eulerian circuit) of a \e directed + ///graph (if there exists) and it converts to the \c Arc type of the digraph. /// - ///For example - ///if the given digraph is Euler (i.e it has only one nontrivial component - ///and the in-degree is equal to the out-degree for all nodes), - ///the following code will put the arcs of \c g - ///to the vector \c et according to an - ///Euler tour of \c g. + ///For example, if the given digraph has an Euler tour (i.e it has only one + ///non-trivial component and the in-degree is equal to the out-degree + ///for all nodes), then the following code will put the arcs of \c g + ///to the vector \c et according to an Euler tour of \c g. ///\code /// std::vector et; - /// for(DiEulerIt e(g),e!=INVALID;++e) + /// for(DiEulerIt e(g); e!=INVALID; ++e) /// et.push_back(e); ///\endcode - ///If \c g is not Euler then the resulted tour will not be full or closed. + ///If \c g has no Euler tour, then the resulted walk will not be closed + ///or not contain all arcs. ///\sa EulerIt template class DiEulerIt @@ -65,53 +63,65 @@ typedef typename GR::InArcIt InArcIt; const GR &g; - typename GR::template NodeMap nedge; + typename GR::template NodeMap narc; std::list euler; public: ///Constructor + ///Constructor. ///\param gr A digraph. - ///\param start The starting point of the tour. If it is not given - /// the tour will start from the first node. + ///\param start The starting point of the tour. If it is not given, + ///the tour will start from the first node that has an outgoing arc. DiEulerIt(const GR &gr, typename GR::Node start = INVALID) - : g(gr), nedge(g) + : g(gr), narc(g) { - if(start==INVALID) start=NodeIt(g); - for(NodeIt n(g);n!=INVALID;++n) nedge[n]=OutArcIt(g,n); - while(nedge[start]!=INVALID) { - euler.push_back(nedge[start]); - Node next=g.target(nedge[start]); - ++nedge[start]; - start=next; + if (start==INVALID) { + NodeIt n(g); + while (n!=INVALID && OutArcIt(g,n)==INVALID) ++n; + start=n; + } + if (start!=INVALID) { + for (NodeIt n(g); n!=INVALID; ++n) narc[n]=OutArcIt(g,n); + while (narc[start]!=INVALID) { + euler.push_back(narc[start]); + Node next=g.target(narc[start]); + ++narc[start]; + start=next; + } } } - ///Arc Conversion + ///Arc conversion operator Arc() { return euler.empty()?INVALID:euler.front(); } + ///Compare with \c INVALID bool operator==(Invalid) { return euler.empty(); } + ///Compare with \c INVALID bool operator!=(Invalid) { return !euler.empty(); } ///Next arc of the tour + + ///Next arc of the tour + /// DiEulerIt &operator++() { Node s=g.target(euler.front()); euler.pop_front(); - //This produces a warning.Strange. - //std::list::iterator next=euler.begin(); typename std::list::iterator next=euler.begin(); - while(nedge[s]!=INVALID) { - euler.insert(next,nedge[s]); - Node n=g.target(nedge[s]); - ++nedge[s]; + while(narc[s]!=INVALID) { + euler.insert(next,narc[s]); + Node n=g.target(narc[s]); + ++narc[s]; s=n; } return *this; } ///Postfix incrementation + /// Postfix incrementation. + /// ///\warning This incrementation - ///returns an \c Arc, not an \ref DiEulerIt, as one may + ///returns an \c Arc, not a \ref DiEulerIt, as one may ///expect. Arc operator++(int) { @@ -121,30 +131,28 @@ } }; - ///Euler iterator for graphs. + ///Euler tour iterator for graphs. /// \ingroup graph_properties - ///This iterator converts to the \c Arc (or \c Edge) - ///type of the digraph and using - ///operator ++, it provides an Euler tour of an undirected - ///digraph (if there exists). + ///This iterator provides an Euler tour (Eulerian circuit) of an + ///\e undirected graph (if there exists) and it converts to the \c Arc + ///and \c Edge types of the graph. /// - ///For example - ///if the given digraph if Euler (i.e it has only one nontrivial component - ///and the degree of each node is even), + ///For example, if the given graph has an Euler tour (i.e it has only one + ///non-trivial component and the degree of each node is even), ///the following code will print the arc IDs according to an ///Euler tour of \c g. ///\code - /// for(EulerIt e(g),e!=INVALID;++e) { + /// for(EulerIt e(g); e!=INVALID; ++e) { /// std::cout << g.id(Edge(e)) << std::eol; /// } ///\endcode - ///Although the iterator provides an Euler tour of an graph, - ///it still returns Arcs in order to indicate the direction of the tour. - ///(But Arc will convert to Edges, of course). + ///Although this iterator is for undirected graphs, it still returns + ///arcs in order to indicate the direction of the tour. + ///(But arcs convert to edges, of course.) /// - ///If \c g is not Euler then the resulted tour will not be full or closed. - ///\sa EulerIt + ///If \c g has no Euler tour, then the resulted walk will not be closed + ///or not contain all edges. template class EulerIt { @@ -157,7 +165,7 @@ typedef typename GR::InArcIt InArcIt; const GR &g; - typename GR::template NodeMap nedge; + typename GR::template NodeMap narc; typename GR::template EdgeMap visited; std::list euler; @@ -165,47 +173,56 @@ ///Constructor - ///\param gr An graph. - ///\param start The starting point of the tour. If it is not given - /// the tour will start from the first node. + ///Constructor. + ///\param gr A graph. + ///\param start The starting point of the tour. If it is not given, + ///the tour will start from the first node that has an incident edge. EulerIt(const GR &gr, typename GR::Node start = INVALID) - : g(gr), nedge(g), visited(g, false) + : g(gr), narc(g), visited(g, false) { - if(start==INVALID) start=NodeIt(g); - for(NodeIt n(g);n!=INVALID;++n) nedge[n]=OutArcIt(g,n); - while(nedge[start]!=INVALID) { - euler.push_back(nedge[start]); - visited[nedge[start]]=true; - Node next=g.target(nedge[start]); - ++nedge[start]; - start=next; - while(nedge[start]!=INVALID && visited[nedge[start]]) ++nedge[start]; + if (start==INVALID) { + NodeIt n(g); + while (n!=INVALID && OutArcIt(g,n)==INVALID) ++n; + start=n; + } + if (start!=INVALID) { + for (NodeIt n(g); n!=INVALID; ++n) narc[n]=OutArcIt(g,n); + while(narc[start]!=INVALID) { + euler.push_back(narc[start]); + visited[narc[start]]=true; + Node next=g.target(narc[start]); + ++narc[start]; + start=next; + while(narc[start]!=INVALID && visited[narc[start]]) ++narc[start]; + } } } - ///Arc Conversion + ///Arc conversion operator Arc() const { return euler.empty()?INVALID:euler.front(); } - ///Arc Conversion + ///Edge conversion operator Edge() const { return euler.empty()?INVALID:euler.front(); } - ///\e + ///Compare with \c INVALID bool operator==(Invalid) const { return euler.empty(); } - ///\e + ///Compare with \c INVALID bool operator!=(Invalid) const { return !euler.empty(); } ///Next arc of the tour + + ///Next arc of the tour + /// EulerIt &operator++() { Node s=g.target(euler.front()); euler.pop_front(); typename std::list::iterator next=euler.begin(); - - while(nedge[s]!=INVALID) { - while(nedge[s]!=INVALID && visited[nedge[s]]) ++nedge[s]; - if(nedge[s]==INVALID) break; + while(narc[s]!=INVALID) { + while(narc[s]!=INVALID && visited[narc[s]]) ++narc[s]; + if(narc[s]==INVALID) break; else { - euler.insert(next,nedge[s]); - visited[nedge[s]]=true; - Node n=g.target(nedge[s]); - ++nedge[s]; + euler.insert(next,narc[s]); + visited[narc[s]]=true; + Node n=g.target(narc[s]); + ++narc[s]; s=n; } } @@ -214,9 +231,10 @@ ///Postfix incrementation - ///\warning This incrementation - ///returns an \c Arc, not an \ref EulerIt, as one may - ///expect. + /// Postfix incrementation. + /// + ///\warning This incrementation returns an \c Arc (which converts to + ///an \c Edge), not an \ref EulerIt, as one may expect. Arc operator++(int) { Arc e=*this; @@ -226,18 +244,23 @@ }; - ///Checks if the graph is Eulerian + ///Check if the given graph is \e Eulerian /// \ingroup graph_properties - ///Checks if the graph is Eulerian. It works for both directed and undirected - ///graphs. - ///\note By definition, a digraph is called \e Eulerian if - ///and only if it is connected and the number of its incoming and outgoing + ///This function checks if the given graph is \e Eulerian. + ///It works for both directed and undirected graphs. + /// + ///By definition, a digraph is called \e Eulerian if + ///and only if it is connected and the number of incoming and outgoing ///arcs are the same for each node. ///Similarly, an undirected graph is called \e Eulerian if - ///and only if it is connected and the number of incident arcs is even - ///for each node. Therefore, there are digraphs which are not Eulerian, - ///but still have an Euler tour. + ///and only if it is connected and the number of incident edges is even + ///for each node. + /// + ///\note There are (di)graphs that are not Eulerian, but still have an + /// Euler tour, since they may contain isolated nodes. + /// + ///\sa DiEulerIt, EulerIt template #ifdef DOXYGEN bool @@ -256,7 +279,7 @@ { for(typename GR::NodeIt n(g);n!=INVALID;++n) if(countInArcs(g,n)!=countOutArcs(g,n)) return false; - return connected(Undirector(g)); + return connected(undirector(g)); } } diff --git a/lemon/matching.h b/lemon/matching.h new file mode 100644 --- /dev/null +++ b/lemon/matching.h @@ -0,0 +1,3244 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2009 + * 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. + * + */ + +#ifndef LEMON_MAX_MATCHING_H +#define LEMON_MAX_MATCHING_H + +#include +#include +#include +#include + +#include +#include +#include +#include + +///\ingroup matching +///\file +///\brief Maximum matching algorithms in general graphs. + +namespace lemon { + + /// \ingroup matching + /// + /// \brief Maximum cardinality matching in general graphs + /// + /// This class implements Edmonds' alternating forest matching algorithm + /// for finding a maximum cardinality matching in a general undirected graph. + /// It can be started from an arbitrary initial matching + /// (the default is the empty one). + /// + /// The dual solution of the problem is a map of the nodes to + /// \ref MaxMatching::Status "Status", having values \c EVEN (or \c D), + /// \c ODD (or \c A) and \c MATCHED (or \c C) defining the Gallai-Edmonds + /// decomposition of the graph. The nodes in \c EVEN/D induce a subgraph + /// 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. The number of the factor-critical components + /// minus the number of barrier nodes is a lower bound on the + /// unmatched nodes, and the matching is optimal if and only if this bound is + /// tight. This decomposition can be obtained using \ref status() or + /// \ref statusMap() after running the algorithm. + /// + /// \tparam GR The undirected graph type the algorithm runs on. + template + class MaxMatching { + public: + + /// The graph type of the algorithm + typedef GR Graph; + /// The type of the matching map + typedef typename Graph::template NodeMap + MatchingMap; + + ///\brief Status constants for Gallai-Edmonds decomposition. + /// + ///These constants are used for indicating the Gallai-Edmonds + ///decomposition of a graph. The nodes with status \c EVEN (or \c D) + ///induce a subgraph with factor-critical components, the nodes with + ///status \c ODD (or \c A) form the canonical barrier, and the nodes + ///with status \c MATCHED (or \c C) induce a subgraph having a + ///perfect matching. + enum Status { + EVEN = 1, ///< = 1. (\c D is an alias for \c EVEN.) + D = 1, + MATCHED = 0, ///< = 0. (\c C is an alias for \c MATCHED.) + C = 0, + ODD = -1, ///< = -1. (\c A is an alias for \c ODD.) + A = -1, + UNMATCHED = -2 ///< = -2. + }; + + /// The type of the status map + typedef typename Graph::template NodeMap StatusMap; + + private: + + TEMPLATE_GRAPH_TYPEDEFS(Graph); + + 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)[node] = arc; + + Node n = node; + while (n != base) { + n = _graph.target((*_matching)[n]); + Arc a = (*_ear)[n]; + n = _graph.target(a); + (*_ear)[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)[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)[_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)[node] = arc; + + Node n = node; + while (n != base) { + n = _graph.target((*_matching)[n]); + Arc a = (*_ear)[n]; + n = _graph.target(a); + (*_ear)[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)[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)[_blossom_set->find(nca)] = nca; + } + + void extendOnArc(const Arc& a) { + Node base = _graph.source(a); + Node odd = _graph.target(a); + + (*_ear)[odd] = _graph.oppositeArc(a); + Node even = _graph.target((*_matching)[odd]); + (*_blossom_rep)[_blossom_set->insert(even)] = even; + (*_status)[odd] = ODD; + (*_status)[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)[odd] = _graph.oppositeArc(a); + (*_status)[odd] = MATCHED; + + Arc arc = (*_matching)[even]; + (*_matching)[even] = a; + + while (arc != INVALID) { + odd = _graph.target(arc); + arc = (*_ear)[odd]; + even = _graph.target(arc); + (*_matching)[odd] = arc; + arc = (*_matching)[even]; + (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]); + } + + for (typename TreeSet::ItemIt it(*_tree_set, tree); + it != INVALID; ++it) { + if ((*_status)[it] == ODD) { + (*_status)[it] = MATCHED; + } else { + int blossom = _blossom_set->find(it); + for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom); + jt != INVALID; ++jt) { + (*_status)[jt] = MATCHED; + } + _blossom_set->eraseClass(blossom); + } + } + _tree_set->eraseClass(tree); + + } + + public: + + /// \brief Constructor + /// + /// 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 + /// \c run() member function.\n + /// If you need better control on the execution, you have to call + /// one of the functions \ref init(), \ref greedyInit() or + /// \ref matchingInit() first, then you can start the algorithm with + /// \ref startSparse() or \ref startDense(). + + ///@{ + + /// \brief Set the initial matching to the empty matching. + /// + /// This function sets the initial matching to the empty matching. + void init() { + createStructures(); + for(NodeIt n(_graph); n != INVALID; ++n) { + (*_matching)[n] = INVALID; + (*_status)[n] = UNMATCHED; + } + } + + /// \brief Find an initial matching in a greedy way. + /// + /// This function finds an initial matching in a greedy way. + void greedyInit() { + createStructures(); + for (NodeIt n(_graph); n != INVALID; ++n) { + (*_matching)[n] = INVALID; + (*_status)[n] = UNMATCHED; + } + 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)[n] = a; + (*_status)[n] = MATCHED; + (*_matching)[v] = _graph.oppositeArc(a); + (*_status)[v] = MATCHED; + break; + } + } + } + } + } + + + /// \brief Initialize the matching from a map. + /// + /// This function initializes the matching from a \c bool valued edge + /// map. This map should have the property that there are no two incident + /// edges with \c true value, i.e. it really contains a matching. + /// \return \c true if the map contains a matching. + template + bool matchingInit(const MatchingMap& matching) { + createStructures(); + + for (NodeIt n(_graph); n != INVALID; ++n) { + (*_matching)[n] = INVALID; + (*_status)[n] = UNMATCHED; + } + for(EdgeIt e(_graph); e!=INVALID; ++e) { + if (matching[e]) { + + Node u = _graph.u(e); + if ((*_matching)[u] != INVALID) return false; + (*_matching)[u] = _graph.direct(e, true); + (*_status)[u] = MATCHED; + + Node v = _graph.v(e); + if ((*_matching)[v] != INVALID) return false; + (*_matching)[v] = _graph.direct(e, false); + (*_status)[v] = MATCHED; + } + } + return true; + } + + /// \brief Start Edmonds' algorithm + /// + /// This function runs the original Edmonds' algorithm. + /// + /// \pre \ref Init(), \ref greedyInit() or \ref matchingInit() must be + /// called before using this function. + 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)[n] = EVEN; + processSparse(n); + } + } + } + + /// \brief Start Edmonds' algorithm with a heuristic improvement + /// for dense graphs + /// + /// This function runs Edmonds' algorithm with a heuristic of postponing + /// shrinks, therefore resulting in a faster algorithm for dense graphs. + /// + /// \pre \ref Init(), \ref greedyInit() or \ref matchingInit() must be + /// called before using this function. + 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)[n] = EVEN; + processDense(n); + } + } + } + + + /// \brief Run Edmonds' algorithm + /// + /// This function runs Edmonds' algorithm. An additional heuristic of + /// postponing shrinks is used for relatively dense graphs + /// (for which m>=2*n holds). + void run() { + if (countEdges(_graph) < 2 * countNodes(_graph)) { + greedyInit(); + startSparse(); + } else { + init(); + startDense(); + } + } + + /// @} + + /// \name Primal Solution + /// Functions to get the primal solution, i.e. the maximum matching. + + /// @{ + + /// \brief Return the size (cardinality) of the matching. + /// + /// This function returns the size (cardinality) of the current matching. + /// After run() it returns the size of the maximum matching in the graph. + int matchingSize() const { + int size = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_matching)[n] != INVALID) { + ++size; + } + } + return size / 2; + } + + /// \brief Return \c true if the given edge is in the matching. + /// + /// This function returns \c true if the given edge is in the current + /// matching. + bool matching(const Edge& edge) const { + return edge == (*_matching)[_graph.u(edge)]; + } + + /// \brief Return the matching arc (or edge) incident to the given node. + /// + /// This function returns the matching arc (or edge) incident to the + /// given node in the current matching or \c INVALID if the node is + /// not covered by the matching. + Arc matching(const Node& n) const { + return (*_matching)[n]; + } + + /// \brief Return a const reference to the matching map. + /// + /// This function returns a const reference to a node map that stores + /// the matching arc (or edge) incident to each node. + const MatchingMap& matchingMap() const { + return *_matching; + } + + /// \brief Return the mate of the given node. + /// + /// This function returns the mate of the given node in the current + /// matching or \c INVALID if the node is not covered by the matching. + Node mate(const Node& n) const { + return (*_matching)[n] != INVALID ? + _graph.target((*_matching)[n]) : INVALID; + } + + /// @} + + /// \name Dual Solution + /// Functions to get the dual solution, i.e. the Gallai-Edmonds + /// decomposition. + + /// @{ + + /// \brief Return the status of the given node in the Edmonds-Gallai + /// decomposition. + /// + /// This function returns the \ref Status "status" of the given node + /// in the Edmonds-Gallai decomposition. + Status status(const Node& n) const { + return (*_status)[n]; + } + + /// \brief Return a const reference to the status map, which stores + /// the Edmonds-Gallai decomposition. + /// + /// This function returns a const reference to a node map that stores the + /// \ref Status "status" of each node in the Edmonds-Gallai decomposition. + const StatusMap& statusMap() const { + return *_status; + } + + /// \brief Return \c true if the given node is in the barrier. + /// + /// This function returns \c true if the given node is in the barrier. + bool barrier(const Node& n) const { + return (*_status)[n] == ODD; + } + + /// @} + + }; + + /// \ingroup matching + /// + /// \brief Weighted matching in general graphs + /// + /// This class provides an efficient implementation of Edmond's + /// maximum weighted 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 a subset of the + /// edges in an undirected graph with maximum overall weight for which + /// each node has at most one incident edge. + /// It can be formulated with the following linear program. + /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f] + /** \f[ \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 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 + /// following. + /** \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] */ + /// + /// The algorithm can be executed with the run() function. + /// After it the matching (the primal solution) and the dual solution + /// can be obtained using the query functions and the + /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, + /// which is able to iterate on the nodes of a blossom. + /// If the value type is integer, then the dual solution is multiplied + /// by \ref MaxWeightedMatching::dualScale "4". + /// + /// \tparam GR The undirected graph type the algorithm runs on. + /// \tparam WM The type edge weight map. The default type is + /// \ref concepts::Graph::EdgeMap "GR::EdgeMap". +#ifdef DOXYGEN + template +#else + template > +#endif + class MaxWeightedMatching { + public: + + /// The graph type of the algorithm + typedef GR Graph; + /// The type of the edge weight map + typedef WM WeightMap; + /// The value type of the edge weights + typedef typename WeightMap::Value Value; + + /// The type of the matching map + typedef typename Graph::template NodeMap + MatchingMap; + + /// \brief Scaling factor for dual solution + /// + /// Scaling factor for dual solution. It is equal to 4 or 1 + /// according to the value type. + static const int dualScale = + std::numeric_limits::is_integer ? 4 : 1; + + private: + + TEMPLATE_GRAPH_TYPEDEFS(Graph); + + typedef typename Graph::template NodeMap NodePotential; + typedef std::vector BlossomNodeList; + + struct BlossomVariable { + int begin, end; + Value value; + + BlossomVariable(int _begin, int _end, Value _value) + : begin(_begin), end(_end), value(_value) {} + + }; + + typedef std::vector BlossomPotential; + + const Graph& _graph; + const WeightMap& _weight; + + MatchingMap* _matching; + + NodePotential* _node_potential; + + BlossomPotential _blossom_potential; + BlossomNodeList _blossom_node_list; + + int _node_num; + int _blossom_num; + + typedef RangeMap IntIntMap; + + enum Status { + EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2 + }; + + typedef HeapUnionFind BlossomSet; + struct BlossomData { + int tree; + Status status; + Arc pred, next; + Value pot, offset; + Node base; + }; + + IntNodeMap *_blossom_index; + BlossomSet *_blossom_set; + RangeMap* _blossom_data; + + IntNodeMap *_node_index; + IntArcMap *_node_heap_index; + + struct NodeData { + + NodeData(IntArcMap& node_heap_index) + : heap(node_heap_index) {} + + int blossom; + Value pot; + BinHeap heap; + std::map heap_index; + + int tree; + }; + + RangeMap* _node_data; + + typedef ExtendFindEnum TreeSet; + + IntIntMap *_tree_set_index; + TreeSet *_tree_set; + + IntNodeMap *_delta1_index; + BinHeap *_delta1; + + IntIntMap *_delta2_index; + BinHeap *_delta2; + + IntEdgeMap *_delta3_index; + BinHeap *_delta3; + + IntIntMap *_delta4_index; + BinHeap *_delta4; + + Value _delta_sum; + + void createStructures() { + _node_num = countNodes(_graph); + _blossom_num = _node_num * 3 / 2; + + if (!_matching) { + _matching = new MatchingMap(_graph); + } + if (!_node_potential) { + _node_potential = new NodePotential(_graph); + } + if (!_blossom_set) { + _blossom_index = new IntNodeMap(_graph); + _blossom_set = new BlossomSet(*_blossom_index); + _blossom_data = new RangeMap(_blossom_num); + } + + if (!_node_index) { + _node_index = new IntNodeMap(_graph); + _node_heap_index = new IntArcMap(_graph); + _node_data = new RangeMap(_node_num, + NodeData(*_node_heap_index)); + } + + if (!_tree_set) { + _tree_set_index = new IntIntMap(_blossom_num); + _tree_set = new TreeSet(*_tree_set_index); + } + if (!_delta1) { + _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 IntEdgeMap(_graph); + _delta3 = new BinHeap(*_delta3_index); + } + if (!_delta4) { + _delta4_index = new IntIntMap(_blossom_num); + _delta4 = new BinHeap(*_delta4_index); + } + } + + void destroyStructures() { + _node_num = countNodes(_graph); + _blossom_num = _node_num * 3 / 2; + + if (_matching) { + delete _matching; + } + if (_node_potential) { + delete _node_potential; + } + if (_blossom_set) { + delete _blossom_index; + delete _blossom_set; + delete _blossom_data; + } + + if (_node_index) { + delete _node_index; + delete _node_heap_index; + delete _node_data; + } + + if (_tree_set) { + delete _tree_set_index; + delete _tree_set; + } + if (_delta1) { + delete _delta1_index; + delete _delta1; + } + if (_delta2) { + delete _delta2_index; + delete _delta2; + } + if (_delta3) { + delete _delta3_index; + delete _delta3; + } + if (_delta4) { + delete _delta4_index; + delete _delta4; + } + } + + void matchedToEven(int blossom, int tree) { + if (_delta2->state(blossom) == _delta2->IN_HEAP) { + _delta2->erase(blossom); + } + + if (!_blossom_set->trivial(blossom)) { + (*_blossom_data)[blossom].pot -= + 2 * (_delta_sum - (*_blossom_data)[blossom].offset); + } + + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); + n != INVALID; ++n) { + + _blossom_set->increase(n, std::numeric_limits::max()); + int ni = (*_node_index)[n]; + + (*_node_data)[ni].heap.clear(); + (*_node_data)[ni].heap_index.clear(); + + (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset; + + _delta1->push(n, (*_node_data)[ni].pot); + + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Node v = _graph.source(e); + int vb = _blossom_set->find(v); + int vi = (*_node_index)[v]; + + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - + dualScale * _weight[e]; + + if ((*_blossom_data)[vb].status == EVEN) { + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { + _delta3->push(e, rw / 2); + } + } else if ((*_blossom_data)[vb].status == UNMATCHED) { + if (_delta3->state(e) != _delta3->IN_HEAP) { + _delta3->push(e, rw); + } + } else { + typename std::map::iterator it = + (*_node_data)[vi].heap_index.find(tree); + + if (it != (*_node_data)[vi].heap_index.end()) { + if ((*_node_data)[vi].heap[it->second] > rw) { + (*_node_data)[vi].heap.replace(it->second, e); + (*_node_data)[vi].heap.decrease(e, rw); + it->second = e; + } + } else { + (*_node_data)[vi].heap.push(e, rw); + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); + } + + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); + + if ((*_blossom_data)[vb].status == MATCHED) { + if (_delta2->state(vb) != _delta2->IN_HEAP) { + _delta2->push(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset){ + _delta2->decrease(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } + } + } + } + } + } + (*_blossom_data)[blossom].offset = 0; + } + + void matchedToOdd(int blossom) { + if (_delta2->state(blossom) == _delta2->IN_HEAP) { + _delta2->erase(blossom); + } + (*_blossom_data)[blossom].offset += _delta_sum; + if (!_blossom_set->trivial(blossom)) { + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + + (*_blossom_data)[blossom].offset); + } + } + + void evenToMatched(int blossom, int tree) { + if (!_blossom_set->trivial(blossom)) { + (*_blossom_data)[blossom].pot += 2 * _delta_sum; + } + + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); + n != INVALID; ++n) { + int ni = (*_node_index)[n]; + (*_node_data)[ni].pot -= _delta_sum; + + _delta1->erase(n); + + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Node v = _graph.source(e); + int vb = _blossom_set->find(v); + int vi = (*_node_index)[v]; + + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - + dualScale * _weight[e]; + + if (vb == blossom) { + if (_delta3->state(e) == _delta3->IN_HEAP) { + _delta3->erase(e); + } + } else if ((*_blossom_data)[vb].status == EVEN) { + + if (_delta3->state(e) == _delta3->IN_HEAP) { + _delta3->erase(e); + } + + int vt = _tree_set->find(vb); + + if (vt != tree) { + + Arc r = _graph.oppositeArc(e); + + typename std::map::iterator it = + (*_node_data)[ni].heap_index.find(vt); + + if (it != (*_node_data)[ni].heap_index.end()) { + if ((*_node_data)[ni].heap[it->second] > rw) { + (*_node_data)[ni].heap.replace(it->second, r); + (*_node_data)[ni].heap.decrease(r, rw); + it->second = r; + } + } else { + (*_node_data)[ni].heap.push(r, rw); + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); + } + + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); + + if (_delta2->state(blossom) != _delta2->IN_HEAP) { + _delta2->push(blossom, _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset); + } else if ((*_delta2)[blossom] > + _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset){ + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset); + } + } + } + + } else if ((*_blossom_data)[vb].status == UNMATCHED) { + if (_delta3->state(e) == _delta3->IN_HEAP) { + _delta3->erase(e); + } + } else { + + typename std::map::iterator it = + (*_node_data)[vi].heap_index.find(tree); + + if (it != (*_node_data)[vi].heap_index.end()) { + (*_node_data)[vi].heap.erase(it->second); + (*_node_data)[vi].heap_index.erase(it); + if ((*_node_data)[vi].heap.empty()) { + _blossom_set->increase(v, std::numeric_limits::max()); + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { + _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); + } + + if ((*_blossom_data)[vb].status == MATCHED) { + if (_blossom_set->classPrio(vb) == + std::numeric_limits::max()) { + _delta2->erase(vb); + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset) { + _delta2->increase(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } + } + } + } + } + } + } + + void oddToMatched(int blossom) { + (*_blossom_data)[blossom].offset -= _delta_sum; + + if (_blossom_set->classPrio(blossom) != + std::numeric_limits::max()) { + _delta2->push(blossom, _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset); + } + + if (!_blossom_set->trivial(blossom)) { + _delta4->erase(blossom); + } + } + + void oddToEven(int blossom, int tree) { + if (!_blossom_set->trivial(blossom)) { + _delta4->erase(blossom); + (*_blossom_data)[blossom].pot -= + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); + } + + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); + n != INVALID; ++n) { + int ni = (*_node_index)[n]; + + _blossom_set->increase(n, std::numeric_limits::max()); + + (*_node_data)[ni].heap.clear(); + (*_node_data)[ni].heap_index.clear(); + (*_node_data)[ni].pot += + 2 * _delta_sum - (*_blossom_data)[blossom].offset; + + _delta1->push(n, (*_node_data)[ni].pot); + + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Node v = _graph.source(e); + int vb = _blossom_set->find(v); + int vi = (*_node_index)[v]; + + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - + dualScale * _weight[e]; + + if ((*_blossom_data)[vb].status == EVEN) { + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { + _delta3->push(e, rw / 2); + } + } else if ((*_blossom_data)[vb].status == UNMATCHED) { + if (_delta3->state(e) != _delta3->IN_HEAP) { + _delta3->push(e, rw); + } + } else { + + typename std::map::iterator it = + (*_node_data)[vi].heap_index.find(tree); + + if (it != (*_node_data)[vi].heap_index.end()) { + if ((*_node_data)[vi].heap[it->second] > rw) { + (*_node_data)[vi].heap.replace(it->second, e); + (*_node_data)[vi].heap.decrease(e, rw); + it->second = e; + } + } else { + (*_node_data)[vi].heap.push(e, rw); + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); + } + + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); + + if ((*_blossom_data)[vb].status == MATCHED) { + if (_delta2->state(vb) != _delta2->IN_HEAP) { + _delta2->push(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset) { + _delta2->decrease(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } + } + } + } + } + } + (*_blossom_data)[blossom].offset = 0; + } + + + void matchedToUnmatched(int blossom) { + if (_delta2->state(blossom) == _delta2->IN_HEAP) { + _delta2->erase(blossom); + } + + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); + n != INVALID; ++n) { + int ni = (*_node_index)[n]; + + _blossom_set->increase(n, std::numeric_limits::max()); + + (*_node_data)[ni].heap.clear(); + (*_node_data)[ni].heap_index.clear(); + + for (OutArcIt e(_graph, n); e != INVALID; ++e) { + Node v = _graph.target(e); + int vb = _blossom_set->find(v); + int vi = (*_node_index)[v]; + + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - + dualScale * _weight[e]; + + if ((*_blossom_data)[vb].status == EVEN) { + if (_delta3->state(e) != _delta3->IN_HEAP) { + _delta3->push(e, rw); + } + } + } + } + } + + void unmatchedToMatched(int blossom) { + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); + n != INVALID; ++n) { + int ni = (*_node_index)[n]; + + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Node v = _graph.source(e); + int vb = _blossom_set->find(v); + int vi = (*_node_index)[v]; + + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - + dualScale * _weight[e]; + + if (vb == blossom) { + if (_delta3->state(e) == _delta3->IN_HEAP) { + _delta3->erase(e); + } + } else if ((*_blossom_data)[vb].status == EVEN) { + + if (_delta3->state(e) == _delta3->IN_HEAP) { + _delta3->erase(e); + } + + int vt = _tree_set->find(vb); + + Arc r = _graph.oppositeArc(e); + + typename std::map::iterator it = + (*_node_data)[ni].heap_index.find(vt); + + if (it != (*_node_data)[ni].heap_index.end()) { + if ((*_node_data)[ni].heap[it->second] > rw) { + (*_node_data)[ni].heap.replace(it->second, r); + (*_node_data)[ni].heap.decrease(r, rw); + it->second = r; + } + } else { + (*_node_data)[ni].heap.push(r, rw); + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); + } + + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); + + if (_delta2->state(blossom) != _delta2->IN_HEAP) { + _delta2->push(blossom, _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset); + } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)- + (*_blossom_data)[blossom].offset){ + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset); + } + } + + } else if ((*_blossom_data)[vb].status == UNMATCHED) { + if (_delta3->state(e) == _delta3->IN_HEAP) { + _delta3->erase(e); + } + } + } + } + } + + void alternatePath(int even, int tree) { + int odd; + + evenToMatched(even, tree); + (*_blossom_data)[even].status = MATCHED; + + while ((*_blossom_data)[even].pred != INVALID) { + odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred)); + (*_blossom_data)[odd].status = MATCHED; + oddToMatched(odd); + (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; + + even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred)); + (*_blossom_data)[even].status = MATCHED; + evenToMatched(even, tree); + (*_blossom_data)[even].next = + _graph.oppositeArc((*_blossom_data)[odd].pred); + } + + } + + void destroyTree(int tree) { + for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { + if ((*_blossom_data)[b].status == EVEN) { + (*_blossom_data)[b].status = MATCHED; + evenToMatched(b, tree); + } else if ((*_blossom_data)[b].status == ODD) { + (*_blossom_data)[b].status = MATCHED; + oddToMatched(b); + } + } + _tree_set->eraseClass(tree); + } + + + void unmatchNode(const Node& node) { + int blossom = _blossom_set->find(node); + int tree = _tree_set->find(blossom); + + alternatePath(blossom, tree); + destroyTree(tree); + + (*_blossom_data)[blossom].status = UNMATCHED; + (*_blossom_data)[blossom].base = node; + matchedToUnmatched(blossom); + } + + + void augmentOnEdge(const Edge& edge) { + + int left = _blossom_set->find(_graph.u(edge)); + int right = _blossom_set->find(_graph.v(edge)); + + if ((*_blossom_data)[left].status == EVEN) { + int left_tree = _tree_set->find(left); + alternatePath(left, left_tree); + destroyTree(left_tree); + } else { + (*_blossom_data)[left].status = MATCHED; + unmatchedToMatched(left); + } + + if ((*_blossom_data)[right].status == EVEN) { + int right_tree = _tree_set->find(right); + alternatePath(right, right_tree); + destroyTree(right_tree); + } else { + (*_blossom_data)[right].status = MATCHED; + unmatchedToMatched(right); + } + + (*_blossom_data)[left].next = _graph.direct(edge, true); + (*_blossom_data)[right].next = _graph.direct(edge, false); + } + + void extendOnArc(const Arc& arc) { + int base = _blossom_set->find(_graph.target(arc)); + int tree = _tree_set->find(base); + + int odd = _blossom_set->find(_graph.source(arc)); + _tree_set->insert(odd, tree); + (*_blossom_data)[odd].status = ODD; + matchedToOdd(odd); + (*_blossom_data)[odd].pred = arc; + + int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next)); + (*_blossom_data)[even].pred = (*_blossom_data)[even].next; + _tree_set->insert(even, tree); + (*_blossom_data)[even].status = EVEN; + matchedToEven(even, tree); + } + + void shrinkOnEdge(const Edge& edge, int tree) { + int nca = -1; + std::vector left_path, right_path; + + { + std::set left_set, right_set; + int left = _blossom_set->find(_graph.u(edge)); + left_path.push_back(left); + left_set.insert(left); + + int right = _blossom_set->find(_graph.v(edge)); + right_path.push_back(right); + right_set.insert(right); + + while (true) { + + if ((*_blossom_data)[left].pred == INVALID) break; + + left = + _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); + left_path.push_back(left); + left = + _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); + left_path.push_back(left); + + left_set.insert(left); + + if (right_set.find(left) != right_set.end()) { + nca = left; + break; + } + + if ((*_blossom_data)[right].pred == INVALID) break; + + right = + _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); + right_path.push_back(right); + right = + _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); + right_path.push_back(right); + + right_set.insert(right); + + if (left_set.find(right) != left_set.end()) { + nca = right; + break; + } + + } + + if (nca == -1) { + if ((*_blossom_data)[left].pred == INVALID) { + nca = right; + while (left_set.find(nca) == left_set.end()) { + nca = + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); + right_path.push_back(nca); + nca = + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); + right_path.push_back(nca); + } + } else { + nca = left; + while (right_set.find(nca) == right_set.end()) { + nca = + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); + left_path.push_back(nca); + nca = + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); + left_path.push_back(nca); + } + } + } + } + + std::vector subblossoms; + Arc prev; + + prev = _graph.direct(edge, true); + for (int i = 0; left_path[i] != nca; i += 2) { + subblossoms.push_back(left_path[i]); + (*_blossom_data)[left_path[i]].next = prev; + _tree_set->erase(left_path[i]); + + subblossoms.push_back(left_path[i + 1]); + (*_blossom_data)[left_path[i + 1]].status = EVEN; + oddToEven(left_path[i + 1], tree); + _tree_set->erase(left_path[i + 1]); + prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred); + } + + int k = 0; + while (right_path[k] != nca) ++k; + + subblossoms.push_back(nca); + (*_blossom_data)[nca].next = prev; + + for (int i = k - 2; i >= 0; i -= 2) { + subblossoms.push_back(right_path[i + 1]); + (*_blossom_data)[right_path[i + 1]].status = EVEN; + oddToEven(right_path[i + 1], tree); + _tree_set->erase(right_path[i + 1]); + + (*_blossom_data)[right_path[i + 1]].next = + (*_blossom_data)[right_path[i + 1]].pred; + + subblossoms.push_back(right_path[i]); + _tree_set->erase(right_path[i]); + } + + int surface = + _blossom_set->join(subblossoms.begin(), subblossoms.end()); + + for (int i = 0; i < int(subblossoms.size()); ++i) { + if (!_blossom_set->trivial(subblossoms[i])) { + (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; + } + (*_blossom_data)[subblossoms[i]].status = MATCHED; + } + + (*_blossom_data)[surface].pot = -2 * _delta_sum; + (*_blossom_data)[surface].offset = 0; + (*_blossom_data)[surface].status = EVEN; + (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; + (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; + + _tree_set->insert(surface, tree); + _tree_set->erase(nca); + } + + void splitBlossom(int blossom) { + Arc next = (*_blossom_data)[blossom].next; + Arc pred = (*_blossom_data)[blossom].pred; + + int tree = _tree_set->find(blossom); + + (*_blossom_data)[blossom].status = MATCHED; + oddToMatched(blossom); + if (_delta2->state(blossom) == _delta2->IN_HEAP) { + _delta2->erase(blossom); + } + + std::vector subblossoms; + _blossom_set->split(blossom, std::back_inserter(subblossoms)); + + Value offset = (*_blossom_data)[blossom].offset; + int b = _blossom_set->find(_graph.source(pred)); + int d = _blossom_set->find(_graph.source(next)); + + int ib = -1, id = -1; + for (int i = 0; i < int(subblossoms.size()); ++i) { + if (subblossoms[i] == b) ib = i; + if (subblossoms[i] == d) id = i; + + (*_blossom_data)[subblossoms[i]].offset = offset; + if (!_blossom_set->trivial(subblossoms[i])) { + (*_blossom_data)[subblossoms[i]].pot -= 2 * offset; + } + if (_blossom_set->classPrio(subblossoms[i]) != + std::numeric_limits::max()) { + _delta2->push(subblossoms[i], + _blossom_set->classPrio(subblossoms[i]) - + (*_blossom_data)[subblossoms[i]].offset); + } + } + + if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) { + for (int i = (id + 1) % subblossoms.size(); + i != ib; i = (i + 2) % subblossoms.size()) { + int sb = subblossoms[i]; + int tb = subblossoms[(i + 1) % subblossoms.size()]; + (*_blossom_data)[sb].next = + _graph.oppositeArc((*_blossom_data)[tb].next); + } + + for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { + int sb = subblossoms[i]; + int tb = subblossoms[(i + 1) % subblossoms.size()]; + int ub = subblossoms[(i + 2) % subblossoms.size()]; + + (*_blossom_data)[sb].status = ODD; + matchedToOdd(sb); + _tree_set->insert(sb, tree); + (*_blossom_data)[sb].pred = pred; + (*_blossom_data)[sb].next = + _graph.oppositeArc((*_blossom_data)[tb].next); + + pred = (*_blossom_data)[ub].next; + + (*_blossom_data)[tb].status = EVEN; + matchedToEven(tb, tree); + _tree_set->insert(tb, tree); + (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; + } + + (*_blossom_data)[subblossoms[id]].status = ODD; + matchedToOdd(subblossoms[id]); + _tree_set->insert(subblossoms[id], tree); + (*_blossom_data)[subblossoms[id]].next = next; + (*_blossom_data)[subblossoms[id]].pred = pred; + + } else { + + for (int i = (ib + 1) % subblossoms.size(); + i != id; i = (i + 2) % subblossoms.size()) { + int sb = subblossoms[i]; + int tb = subblossoms[(i + 1) % subblossoms.size()]; + (*_blossom_data)[sb].next = + _graph.oppositeArc((*_blossom_data)[tb].next); + } + + for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { + int sb = subblossoms[i]; + int tb = subblossoms[(i + 1) % subblossoms.size()]; + int ub = subblossoms[(i + 2) % subblossoms.size()]; + + (*_blossom_data)[sb].status = ODD; + matchedToOdd(sb); + _tree_set->insert(sb, tree); + (*_blossom_data)[sb].next = next; + (*_blossom_data)[sb].pred = + _graph.oppositeArc((*_blossom_data)[tb].next); + + (*_blossom_data)[tb].status = EVEN; + matchedToEven(tb, tree); + _tree_set->insert(tb, tree); + (*_blossom_data)[tb].pred = + (*_blossom_data)[tb].next = + _graph.oppositeArc((*_blossom_data)[ub].next); + next = (*_blossom_data)[ub].next; + } + + (*_blossom_data)[subblossoms[ib]].status = ODD; + matchedToOdd(subblossoms[ib]); + _tree_set->insert(subblossoms[ib], tree); + (*_blossom_data)[subblossoms[ib]].next = next; + (*_blossom_data)[subblossoms[ib]].pred = pred; + } + _tree_set->erase(blossom); + } + + void extractBlossom(int blossom, const Node& base, const Arc& matching) { + if (_blossom_set->trivial(blossom)) { + int bi = (*_node_index)[base]; + Value pot = (*_node_data)[bi].pot; + + (*_matching)[base] = matching; + _blossom_node_list.push_back(base); + (*_node_potential)[base] = pot; + } else { + + Value pot = (*_blossom_data)[blossom].pot; + int bn = _blossom_node_list.size(); + + std::vector subblossoms; + _blossom_set->split(blossom, std::back_inserter(subblossoms)); + int b = _blossom_set->find(base); + int ib = -1; + for (int i = 0; i < int(subblossoms.size()); ++i) { + if (subblossoms[i] == b) { ib = i; break; } + } + + for (int i = 1; i < int(subblossoms.size()); i += 2) { + int sb = subblossoms[(ib + i) % subblossoms.size()]; + int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; + + Arc m = (*_blossom_data)[tb].next; + extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m)); + extractBlossom(tb, _graph.source(m), m); + } + extractBlossom(subblossoms[ib], base, matching); + + int en = _blossom_node_list.size(); + + _blossom_potential.push_back(BlossomVariable(bn, en, pot)); + } + } + + void extractMatching() { + std::vector blossoms; + for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { + blossoms.push_back(c); + } + + for (int i = 0; i < int(blossoms.size()); ++i) { + if ((*_blossom_data)[blossoms[i]].status == MATCHED) { + + Value offset = (*_blossom_data)[blossoms[i]].offset; + (*_blossom_data)[blossoms[i]].pot += 2 * offset; + for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); + n != INVALID; ++n) { + (*_node_data)[(*_node_index)[n]].pot -= offset; + } + + Arc matching = (*_blossom_data)[blossoms[i]].next; + Node base = _graph.source(matching); + extractBlossom(blossoms[i], base, matching); + } else { + Node base = (*_blossom_data)[blossoms[i]].base; + extractBlossom(blossoms[i], base, INVALID); + } + } + } + + public: + + /// \brief Constructor + /// + /// Constructor. + MaxWeightedMatching(const Graph& graph, const WeightMap& weight) + : _graph(graph), _weight(weight), _matching(0), + _node_potential(0), _blossom_potential(), _blossom_node_list(), + _node_num(0), _blossom_num(0), + + _blossom_index(0), _blossom_set(0), _blossom_data(0), + _node_index(0), _node_heap_index(0), _node_data(0), + _tree_set_index(0), _tree_set(0), + + _delta1_index(0), _delta1(0), + _delta2_index(0), _delta2(0), + _delta3_index(0), _delta3(0), + _delta4_index(0), _delta4(0), + + _delta_sum() {} + + ~MaxWeightedMatching() { + destroyStructures(); + } + + /// \name Execution Control + /// The simplest way to execute the algorithm is to use the + /// \ref run() member function. + + ///@{ + + /// \brief Initialize the algorithm + /// + /// This function initializes the algorithm. + void init() { + createStructures(); + + for (ArcIt e(_graph); e != INVALID; ++e) { + (*_node_heap_index)[e] = BinHeap::PRE_HEAP; + } + for (NodeIt n(_graph); n != INVALID; ++n) { + (*_delta1_index)[n] = _delta1->PRE_HEAP; + } + for (EdgeIt e(_graph); e != INVALID; ++e) { + (*_delta3_index)[e] = _delta3->PRE_HEAP; + } + for (int i = 0; i < _blossom_num; ++i) { + (*_delta2_index)[i] = _delta2->PRE_HEAP; + (*_delta4_index)[i] = _delta4->PRE_HEAP; + } + + int index = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + Value max = 0; + for (OutArcIt e(_graph, n); e != INVALID; ++e) { + if (_graph.target(e) == n) continue; + if ((dualScale * _weight[e]) / 2 > max) { + max = (dualScale * _weight[e]) / 2; + } + } + (*_node_index)[n] = index; + (*_node_data)[index].pot = max; + _delta1->push(n, max); + int blossom = + _blossom_set->insert(n, std::numeric_limits::max()); + + _tree_set->insert(blossom); + + (*_blossom_data)[blossom].status = EVEN; + (*_blossom_data)[blossom].pred = INVALID; + (*_blossom_data)[blossom].next = INVALID; + (*_blossom_data)[blossom].pot = 0; + (*_blossom_data)[blossom].offset = 0; + ++index; + } + for (EdgeIt e(_graph); e != INVALID; ++e) { + int si = (*_node_index)[_graph.u(e)]; + int ti = (*_node_index)[_graph.v(e)]; + if (_graph.u(e) != _graph.v(e)) { + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - + dualScale * _weight[e]) / 2); + } + } + } + + /// \brief Start the algorithm + /// + /// This function starts the algorithm. + /// + /// \pre \ref init() must be called before using this function. + void start() { + enum OpType { + D1, D2, D3, D4 + }; + + int unmatched = _node_num; + while (unmatched > 0) { + Value d1 = !_delta1->empty() ? + _delta1->prio() : std::numeric_limits::max(); + + Value d2 = !_delta2->empty() ? + _delta2->prio() : std::numeric_limits::max(); + + Value d3 = !_delta3->empty() ? + _delta3->prio() : std::numeric_limits::max(); + + Value d4 = !_delta4->empty() ? + _delta4->prio() : std::numeric_limits::max(); + + _delta_sum = d1; OpType ot = D1; + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } + if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } + if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } + + + switch (ot) { + case D1: + { + Node n = _delta1->top(); + unmatchNode(n); + --unmatched; + } + break; + case D2: + { + int blossom = _delta2->top(); + Node n = _blossom_set->classTop(blossom); + Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); + extendOnArc(e); + } + break; + case D3: + { + Edge e = _delta3->top(); + + int left_blossom = _blossom_set->find(_graph.u(e)); + int right_blossom = _blossom_set->find(_graph.v(e)); + + if (left_blossom == right_blossom) { + _delta3->pop(); + } else { + int left_tree; + if ((*_blossom_data)[left_blossom].status == EVEN) { + left_tree = _tree_set->find(left_blossom); + } else { + left_tree = -1; + ++unmatched; + } + int right_tree; + if ((*_blossom_data)[right_blossom].status == EVEN) { + right_tree = _tree_set->find(right_blossom); + } else { + right_tree = -1; + ++unmatched; + } + + if (left_tree == right_tree) { + shrinkOnEdge(e, left_tree); + } else { + augmentOnEdge(e); + unmatched -= 2; + } + } + } break; + case D4: + splitBlossom(_delta4->top()); + break; + } + } + extractMatching(); + } + + /// \brief Run the algorithm. + /// + /// This method runs the \c %MaxWeightedMatching algorithm. + /// + /// \note mwm.run() is just a shortcut of the following code. + /// \code + /// mwm.init(); + /// mwm.start(); + /// \endcode + void run() { + init(); + start(); + } + + /// @} + + /// \name Primal Solution + /// Functions to get the primal solution, i.e. the maximum weighted + /// matching.\n + /// Either \ref run() or \ref start() function should be called before + /// using them. + + /// @{ + + /// \brief Return the weight of the matching. + /// + /// This function returns the weight of the found matching. + /// + /// \pre Either run() or start() must be called before using this function. + Value matchingWeight() const { + Value sum = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_matching)[n] != INVALID) { + sum += _weight[(*_matching)[n]]; + } + } + return sum /= 2; + } + + /// \brief Return the size (cardinality) of the matching. + /// + /// This function returns the size (cardinality) of the found matching. + /// + /// \pre Either run() or start() must be called before using this function. + int matchingSize() const { + int num = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_matching)[n] != INVALID) { + ++num; + } + } + return num /= 2; + } + + /// \brief Return \c true if the given edge is in the matching. + /// + /// This function returns \c true if the given edge is in the found + /// matching. + /// + /// \pre Either run() or start() must be called before using this function. + bool matching(const Edge& edge) const { + return edge == (*_matching)[_graph.u(edge)]; + } + + /// \brief Return the matching arc (or edge) incident to the given node. + /// + /// This function returns the matching arc (or edge) incident to the + /// given node in the found matching or \c INVALID if the node is + /// not covered by the matching. + /// + /// \pre Either run() or start() must be called before using this function. + Arc matching(const Node& node) const { + return (*_matching)[node]; + } + + /// \brief Return a const reference to the matching map. + /// + /// This function returns a const reference to a node map that stores + /// the matching arc (or edge) incident to each node. + const MatchingMap& matchingMap() const { + return *_matching; + } + + /// \brief Return the mate of the given node. + /// + /// This function returns the mate of the given node in the found + /// matching or \c INVALID if the node is not covered by the matching. + /// + /// \pre Either run() or start() must be called before using this function. + Node mate(const Node& node) const { + return (*_matching)[node] != INVALID ? + _graph.target((*_matching)[node]) : INVALID; + } + + /// @} + + /// \name Dual Solution + /// Functions to get the dual solution.\n + /// Either \ref run() or \ref start() function should be called before + /// using them. + + /// @{ + + /// \brief Return the value of the dual solution. + /// + /// This function returns the value of the dual solution. + /// It should be equal to the primal value scaled by \ref dualScale + /// "dual scale". + /// + /// \pre Either run() or start() must be called before using this function. + Value dualValue() const { + Value sum = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + sum += nodeValue(n); + } + for (int i = 0; i < blossomNum(); ++i) { + sum += blossomValue(i) * (blossomSize(i) / 2); + } + return sum; + } + + /// \brief Return the dual value (potential) of the given node. + /// + /// This function returns the dual value (potential) of the given node. + /// + /// \pre Either run() or start() must be called before using this function. + Value nodeValue(const Node& n) const { + return (*_node_potential)[n]; + } + + /// \brief Return the number of the blossoms in the basis. + /// + /// This function returns the number of the blossoms in the basis. + /// + /// \pre Either run() or start() must be called before using this function. + /// \see BlossomIt + int blossomNum() const { + return _blossom_potential.size(); + } + + /// \brief Return the number of the nodes in the given blossom. + /// + /// This function returns the number of the nodes in the given blossom. + /// + /// \pre Either run() or start() must be called before using this function. + /// \see BlossomIt + int blossomSize(int k) const { + return _blossom_potential[k].end - _blossom_potential[k].begin; + } + + /// \brief Return the dual value (ptential) of the given blossom. + /// + /// This function returns the dual value (ptential) of the given blossom. + /// + /// \pre Either run() or start() must be called before using this function. + Value blossomValue(int k) const { + return _blossom_potential[k].value; + } + + /// \brief Iterator for obtaining the nodes of a blossom. + /// + /// This class provides an iterator for obtaining the nodes of the + /// given blossom. It lists a subset of the nodes. + /// Before using this iterator, you must allocate a + /// MaxWeightedMatching class and execute it. + class BlossomIt { + public: + + /// \brief Constructor. + /// + /// Constructor to get the nodes of the given variable. + /// + /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or + /// \ref MaxWeightedMatching::start() "algorithm.start()" must be + /// called before initializing this iterator. + BlossomIt(const MaxWeightedMatching& algorithm, int variable) + : _algorithm(&algorithm) + { + _index = _algorithm->_blossom_potential[variable].begin; + _last = _algorithm->_blossom_potential[variable].end; + } + + /// \brief Conversion to \c Node. + /// + /// Conversion to \c Node. + operator Node() const { + return _algorithm->_blossom_node_list[_index]; + } + + /// \brief Increment operator. + /// + /// Increment operator. + BlossomIt& operator++() { + ++_index; + return *this; + } + + /// \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; + int _last; + int _index; + }; + + /// @} + + }; + + /// \ingroup matching + /// + /// \brief Weighted perfect matching in general graphs + /// + /// This class provides an efficient implementation of Edmond's + /// 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 perfect matching problem is to find a subset of + /// the edges in an undirected graph with maximum overall weight for which + /// each node has exactly one incident edge. + /// It can be formulated with the following linear program. + /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f] + /** \f[ \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 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 + /// following. + /** \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] */ + /// + /// The algorithm can be executed with the run() function. + /// After it the matching (the primal solution) and the dual solution + /// can be obtained using the query functions and the + /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, + /// which is able to iterate on the nodes of a blossom. + /// If the value type is integer, then the dual solution is multiplied + /// by \ref MaxWeightedMatching::dualScale "4". + /// + /// \tparam GR The undirected graph type the algorithm runs on. + /// \tparam WM The type edge weight map. The default type is + /// \ref concepts::Graph::EdgeMap "GR::EdgeMap". +#ifdef DOXYGEN + template +#else + template > +#endif + class MaxWeightedPerfectMatching { + public: + + /// The graph type of the algorithm + typedef GR Graph; + /// The type of the edge weight map + typedef WM WeightMap; + /// The value type of the edge weights + typedef typename WeightMap::Value Value; + + /// \brief Scaling factor for dual solution + /// + /// Scaling factor for dual solution, it is equal to 4 or 1 + /// according to the value type. + static const int dualScale = + std::numeric_limits::is_integer ? 4 : 1; + + /// The type of the matching map + typedef typename Graph::template NodeMap + MatchingMap; + + private: + + TEMPLATE_GRAPH_TYPEDEFS(Graph); + + typedef typename Graph::template NodeMap NodePotential; + typedef std::vector BlossomNodeList; + + struct BlossomVariable { + int begin, end; + Value value; + + BlossomVariable(int _begin, int _end, Value _value) + : begin(_begin), end(_end), value(_value) {} + + }; + + typedef std::vector BlossomPotential; + + const Graph& _graph; + const WeightMap& _weight; + + MatchingMap* _matching; + + NodePotential* _node_potential; + + BlossomPotential _blossom_potential; + BlossomNodeList _blossom_node_list; + + int _node_num; + int _blossom_num; + + typedef RangeMap IntIntMap; + + enum Status { + EVEN = -1, MATCHED = 0, ODD = 1 + }; + + typedef HeapUnionFind BlossomSet; + struct BlossomData { + int tree; + Status status; + Arc pred, next; + Value pot, offset; + }; + + IntNodeMap *_blossom_index; + BlossomSet *_blossom_set; + RangeMap* _blossom_data; + + IntNodeMap *_node_index; + IntArcMap *_node_heap_index; + + struct NodeData { + + NodeData(IntArcMap& node_heap_index) + : heap(node_heap_index) {} + + int blossom; + Value pot; + BinHeap heap; + std::map heap_index; + + int tree; + }; + + RangeMap* _node_data; + + typedef ExtendFindEnum TreeSet; + + IntIntMap *_tree_set_index; + TreeSet *_tree_set; + + IntIntMap *_delta2_index; + BinHeap *_delta2; + + IntEdgeMap *_delta3_index; + BinHeap *_delta3; + + IntIntMap *_delta4_index; + BinHeap *_delta4; + + Value _delta_sum; + + void createStructures() { + _node_num = countNodes(_graph); + _blossom_num = _node_num * 3 / 2; + + if (!_matching) { + _matching = new MatchingMap(_graph); + } + if (!_node_potential) { + _node_potential = new NodePotential(_graph); + } + if (!_blossom_set) { + _blossom_index = new IntNodeMap(_graph); + _blossom_set = new BlossomSet(*_blossom_index); + _blossom_data = new RangeMap(_blossom_num); + } + + if (!_node_index) { + _node_index = new IntNodeMap(_graph); + _node_heap_index = new IntArcMap(_graph); + _node_data = new RangeMap(_node_num, + NodeData(*_node_heap_index)); + } + + if (!_tree_set) { + _tree_set_index = new IntIntMap(_blossom_num); + _tree_set = new TreeSet(*_tree_set_index); + } + if (!_delta2) { + _delta2_index = new IntIntMap(_blossom_num); + _delta2 = new BinHeap(*_delta2_index); + } + if (!_delta3) { + _delta3_index = new IntEdgeMap(_graph); + _delta3 = new BinHeap(*_delta3_index); + } + if (!_delta4) { + _delta4_index = new IntIntMap(_blossom_num); + _delta4 = new BinHeap(*_delta4_index); + } + } + + void destroyStructures() { + _node_num = countNodes(_graph); + _blossom_num = _node_num * 3 / 2; + + if (_matching) { + delete _matching; + } + if (_node_potential) { + delete _node_potential; + } + if (_blossom_set) { + delete _blossom_index; + delete _blossom_set; + delete _blossom_data; + } + + if (_node_index) { + delete _node_index; + delete _node_heap_index; + delete _node_data; + } + + if (_tree_set) { + delete _tree_set_index; + delete _tree_set; + } + if (_delta2) { + delete _delta2_index; + delete _delta2; + } + if (_delta3) { + delete _delta3_index; + delete _delta3; + } + if (_delta4) { + delete _delta4_index; + delete _delta4; + } + } + + void matchedToEven(int blossom, int tree) { + if (_delta2->state(blossom) == _delta2->IN_HEAP) { + _delta2->erase(blossom); + } + + if (!_blossom_set->trivial(blossom)) { + (*_blossom_data)[blossom].pot -= + 2 * (_delta_sum - (*_blossom_data)[blossom].offset); + } + + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); + n != INVALID; ++n) { + + _blossom_set->increase(n, std::numeric_limits::max()); + int ni = (*_node_index)[n]; + + (*_node_data)[ni].heap.clear(); + (*_node_data)[ni].heap_index.clear(); + + (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset; + + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Node v = _graph.source(e); + int vb = _blossom_set->find(v); + int vi = (*_node_index)[v]; + + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - + dualScale * _weight[e]; + + if ((*_blossom_data)[vb].status == EVEN) { + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { + _delta3->push(e, rw / 2); + } + } else { + typename std::map::iterator it = + (*_node_data)[vi].heap_index.find(tree); + + if (it != (*_node_data)[vi].heap_index.end()) { + if ((*_node_data)[vi].heap[it->second] > rw) { + (*_node_data)[vi].heap.replace(it->second, e); + (*_node_data)[vi].heap.decrease(e, rw); + it->second = e; + } + } else { + (*_node_data)[vi].heap.push(e, rw); + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); + } + + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); + + if ((*_blossom_data)[vb].status == MATCHED) { + if (_delta2->state(vb) != _delta2->IN_HEAP) { + _delta2->push(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset){ + _delta2->decrease(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } + } + } + } + } + } + (*_blossom_data)[blossom].offset = 0; + } + + void matchedToOdd(int blossom) { + if (_delta2->state(blossom) == _delta2->IN_HEAP) { + _delta2->erase(blossom); + } + (*_blossom_data)[blossom].offset += _delta_sum; + if (!_blossom_set->trivial(blossom)) { + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + + (*_blossom_data)[blossom].offset); + } + } + + void evenToMatched(int blossom, int tree) { + if (!_blossom_set->trivial(blossom)) { + (*_blossom_data)[blossom].pot += 2 * _delta_sum; + } + + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); + n != INVALID; ++n) { + int ni = (*_node_index)[n]; + (*_node_data)[ni].pot -= _delta_sum; + + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Node v = _graph.source(e); + int vb = _blossom_set->find(v); + int vi = (*_node_index)[v]; + + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - + dualScale * _weight[e]; + + if (vb == blossom) { + if (_delta3->state(e) == _delta3->IN_HEAP) { + _delta3->erase(e); + } + } else if ((*_blossom_data)[vb].status == EVEN) { + + if (_delta3->state(e) == _delta3->IN_HEAP) { + _delta3->erase(e); + } + + int vt = _tree_set->find(vb); + + if (vt != tree) { + + Arc r = _graph.oppositeArc(e); + + typename std::map::iterator it = + (*_node_data)[ni].heap_index.find(vt); + + if (it != (*_node_data)[ni].heap_index.end()) { + if ((*_node_data)[ni].heap[it->second] > rw) { + (*_node_data)[ni].heap.replace(it->second, r); + (*_node_data)[ni].heap.decrease(r, rw); + it->second = r; + } + } else { + (*_node_data)[ni].heap.push(r, rw); + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); + } + + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); + + if (_delta2->state(blossom) != _delta2->IN_HEAP) { + _delta2->push(blossom, _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset); + } else if ((*_delta2)[blossom] > + _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset){ + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset); + } + } + } + } else { + + typename std::map::iterator it = + (*_node_data)[vi].heap_index.find(tree); + + if (it != (*_node_data)[vi].heap_index.end()) { + (*_node_data)[vi].heap.erase(it->second); + (*_node_data)[vi].heap_index.erase(it); + if ((*_node_data)[vi].heap.empty()) { + _blossom_set->increase(v, std::numeric_limits::max()); + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { + _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); + } + + if ((*_blossom_data)[vb].status == MATCHED) { + if (_blossom_set->classPrio(vb) == + std::numeric_limits::max()) { + _delta2->erase(vb); + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset) { + _delta2->increase(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } + } + } + } + } + } + } + + void oddToMatched(int blossom) { + (*_blossom_data)[blossom].offset -= _delta_sum; + + if (_blossom_set->classPrio(blossom) != + std::numeric_limits::max()) { + _delta2->push(blossom, _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset); + } + + if (!_blossom_set->trivial(blossom)) { + _delta4->erase(blossom); + } + } + + void oddToEven(int blossom, int tree) { + if (!_blossom_set->trivial(blossom)) { + _delta4->erase(blossom); + (*_blossom_data)[blossom].pot -= + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); + } + + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); + n != INVALID; ++n) { + int ni = (*_node_index)[n]; + + _blossom_set->increase(n, std::numeric_limits::max()); + + (*_node_data)[ni].heap.clear(); + (*_node_data)[ni].heap_index.clear(); + (*_node_data)[ni].pot += + 2 * _delta_sum - (*_blossom_data)[blossom].offset; + + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Node v = _graph.source(e); + int vb = _blossom_set->find(v); + int vi = (*_node_index)[v]; + + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - + dualScale * _weight[e]; + + if ((*_blossom_data)[vb].status == EVEN) { + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { + _delta3->push(e, rw / 2); + } + } else { + + typename std::map::iterator it = + (*_node_data)[vi].heap_index.find(tree); + + if (it != (*_node_data)[vi].heap_index.end()) { + if ((*_node_data)[vi].heap[it->second] > rw) { + (*_node_data)[vi].heap.replace(it->second, e); + (*_node_data)[vi].heap.decrease(e, rw); + it->second = e; + } + } else { + (*_node_data)[vi].heap.push(e, rw); + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); + } + + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); + + if ((*_blossom_data)[vb].status == MATCHED) { + if (_delta2->state(vb) != _delta2->IN_HEAP) { + _delta2->push(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset) { + _delta2->decrease(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } + } + } + } + } + } + (*_blossom_data)[blossom].offset = 0; + } + + void alternatePath(int even, int tree) { + int odd; + + evenToMatched(even, tree); + (*_blossom_data)[even].status = MATCHED; + + while ((*_blossom_data)[even].pred != INVALID) { + odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred)); + (*_blossom_data)[odd].status = MATCHED; + oddToMatched(odd); + (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; + + even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred)); + (*_blossom_data)[even].status = MATCHED; + evenToMatched(even, tree); + (*_blossom_data)[even].next = + _graph.oppositeArc((*_blossom_data)[odd].pred); + } + + } + + void destroyTree(int tree) { + for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { + if ((*_blossom_data)[b].status == EVEN) { + (*_blossom_data)[b].status = MATCHED; + evenToMatched(b, tree); + } else if ((*_blossom_data)[b].status == ODD) { + (*_blossom_data)[b].status = MATCHED; + oddToMatched(b); + } + } + _tree_set->eraseClass(tree); + } + + 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); + destroyTree(left_tree); + + int right_tree = _tree_set->find(right); + alternatePath(right, right_tree); + destroyTree(right_tree); + + (*_blossom_data)[left].next = _graph.direct(edge, true); + (*_blossom_data)[right].next = _graph.direct(edge, false); + } + + void extendOnArc(const Arc& arc) { + int base = _blossom_set->find(_graph.target(arc)); + int tree = _tree_set->find(base); + + int odd = _blossom_set->find(_graph.source(arc)); + _tree_set->insert(odd, tree); + (*_blossom_data)[odd].status = ODD; + matchedToOdd(odd); + (*_blossom_data)[odd].pred = arc; + + int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next)); + (*_blossom_data)[even].pred = (*_blossom_data)[even].next; + _tree_set->insert(even, tree); + (*_blossom_data)[even].status = EVEN; + matchedToEven(even, tree); + } + + void shrinkOnEdge(const Edge& edge, int tree) { + int nca = -1; + std::vector left_path, right_path; + + { + std::set left_set, right_set; + int left = _blossom_set->find(_graph.u(edge)); + left_path.push_back(left); + left_set.insert(left); + + int right = _blossom_set->find(_graph.v(edge)); + right_path.push_back(right); + right_set.insert(right); + + while (true) { + + if ((*_blossom_data)[left].pred == INVALID) break; + + left = + _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); + left_path.push_back(left); + left = + _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); + left_path.push_back(left); + + left_set.insert(left); + + if (right_set.find(left) != right_set.end()) { + nca = left; + break; + } + + if ((*_blossom_data)[right].pred == INVALID) break; + + right = + _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); + right_path.push_back(right); + right = + _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); + right_path.push_back(right); + + right_set.insert(right); + + if (left_set.find(right) != left_set.end()) { + nca = right; + break; + } + + } + + if (nca == -1) { + if ((*_blossom_data)[left].pred == INVALID) { + nca = right; + while (left_set.find(nca) == left_set.end()) { + nca = + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); + right_path.push_back(nca); + nca = + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); + right_path.push_back(nca); + } + } else { + nca = left; + while (right_set.find(nca) == right_set.end()) { + nca = + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); + left_path.push_back(nca); + nca = + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); + left_path.push_back(nca); + } + } + } + } + + std::vector subblossoms; + Arc prev; + + prev = _graph.direct(edge, true); + for (int i = 0; left_path[i] != nca; i += 2) { + subblossoms.push_back(left_path[i]); + (*_blossom_data)[left_path[i]].next = prev; + _tree_set->erase(left_path[i]); + + subblossoms.push_back(left_path[i + 1]); + (*_blossom_data)[left_path[i + 1]].status = EVEN; + oddToEven(left_path[i + 1], tree); + _tree_set->erase(left_path[i + 1]); + prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred); + } + + int k = 0; + while (right_path[k] != nca) ++k; + + subblossoms.push_back(nca); + (*_blossom_data)[nca].next = prev; + + for (int i = k - 2; i >= 0; i -= 2) { + subblossoms.push_back(right_path[i + 1]); + (*_blossom_data)[right_path[i + 1]].status = EVEN; + oddToEven(right_path[i + 1], tree); + _tree_set->erase(right_path[i + 1]); + + (*_blossom_data)[right_path[i + 1]].next = + (*_blossom_data)[right_path[i + 1]].pred; + + subblossoms.push_back(right_path[i]); + _tree_set->erase(right_path[i]); + } + + int surface = + _blossom_set->join(subblossoms.begin(), subblossoms.end()); + + for (int i = 0; i < int(subblossoms.size()); ++i) { + if (!_blossom_set->trivial(subblossoms[i])) { + (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; + } + (*_blossom_data)[subblossoms[i]].status = MATCHED; + } + + (*_blossom_data)[surface].pot = -2 * _delta_sum; + (*_blossom_data)[surface].offset = 0; + (*_blossom_data)[surface].status = EVEN; + (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; + (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; + + _tree_set->insert(surface, tree); + _tree_set->erase(nca); + } + + void splitBlossom(int blossom) { + Arc next = (*_blossom_data)[blossom].next; + Arc pred = (*_blossom_data)[blossom].pred; + + int tree = _tree_set->find(blossom); + + (*_blossom_data)[blossom].status = MATCHED; + oddToMatched(blossom); + if (_delta2->state(blossom) == _delta2->IN_HEAP) { + _delta2->erase(blossom); + } + + std::vector subblossoms; + _blossom_set->split(blossom, std::back_inserter(subblossoms)); + + Value offset = (*_blossom_data)[blossom].offset; + int b = _blossom_set->find(_graph.source(pred)); + int d = _blossom_set->find(_graph.source(next)); + + int ib = -1, id = -1; + for (int i = 0; i < int(subblossoms.size()); ++i) { + if (subblossoms[i] == b) ib = i; + if (subblossoms[i] == d) id = i; + + (*_blossom_data)[subblossoms[i]].offset = offset; + if (!_blossom_set->trivial(subblossoms[i])) { + (*_blossom_data)[subblossoms[i]].pot -= 2 * offset; + } + if (_blossom_set->classPrio(subblossoms[i]) != + std::numeric_limits::max()) { + _delta2->push(subblossoms[i], + _blossom_set->classPrio(subblossoms[i]) - + (*_blossom_data)[subblossoms[i]].offset); + } + } + + if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) { + for (int i = (id + 1) % subblossoms.size(); + i != ib; i = (i + 2) % subblossoms.size()) { + int sb = subblossoms[i]; + int tb = subblossoms[(i + 1) % subblossoms.size()]; + (*_blossom_data)[sb].next = + _graph.oppositeArc((*_blossom_data)[tb].next); + } + + for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { + int sb = subblossoms[i]; + int tb = subblossoms[(i + 1) % subblossoms.size()]; + int ub = subblossoms[(i + 2) % subblossoms.size()]; + + (*_blossom_data)[sb].status = ODD; + matchedToOdd(sb); + _tree_set->insert(sb, tree); + (*_blossom_data)[sb].pred = pred; + (*_blossom_data)[sb].next = + _graph.oppositeArc((*_blossom_data)[tb].next); + + pred = (*_blossom_data)[ub].next; + + (*_blossom_data)[tb].status = EVEN; + matchedToEven(tb, tree); + _tree_set->insert(tb, tree); + (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; + } + + (*_blossom_data)[subblossoms[id]].status = ODD; + matchedToOdd(subblossoms[id]); + _tree_set->insert(subblossoms[id], tree); + (*_blossom_data)[subblossoms[id]].next = next; + (*_blossom_data)[subblossoms[id]].pred = pred; + + } else { + + for (int i = (ib + 1) % subblossoms.size(); + i != id; i = (i + 2) % subblossoms.size()) { + int sb = subblossoms[i]; + int tb = subblossoms[(i + 1) % subblossoms.size()]; + (*_blossom_data)[sb].next = + _graph.oppositeArc((*_blossom_data)[tb].next); + } + + for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { + int sb = subblossoms[i]; + int tb = subblossoms[(i + 1) % subblossoms.size()]; + int ub = subblossoms[(i + 2) % subblossoms.size()]; + + (*_blossom_data)[sb].status = ODD; + matchedToOdd(sb); + _tree_set->insert(sb, tree); + (*_blossom_data)[sb].next = next; + (*_blossom_data)[sb].pred = + _graph.oppositeArc((*_blossom_data)[tb].next); + + (*_blossom_data)[tb].status = EVEN; + matchedToEven(tb, tree); + _tree_set->insert(tb, tree); + (*_blossom_data)[tb].pred = + (*_blossom_data)[tb].next = + _graph.oppositeArc((*_blossom_data)[ub].next); + next = (*_blossom_data)[ub].next; + } + + (*_blossom_data)[subblossoms[ib]].status = ODD; + matchedToOdd(subblossoms[ib]); + _tree_set->insert(subblossoms[ib], tree); + (*_blossom_data)[subblossoms[ib]].next = next; + (*_blossom_data)[subblossoms[ib]].pred = pred; + } + _tree_set->erase(blossom); + } + + void extractBlossom(int blossom, const Node& base, const Arc& matching) { + if (_blossom_set->trivial(blossom)) { + int bi = (*_node_index)[base]; + Value pot = (*_node_data)[bi].pot; + + (*_matching)[base] = matching; + _blossom_node_list.push_back(base); + (*_node_potential)[base] = pot; + } else { + + Value pot = (*_blossom_data)[blossom].pot; + int bn = _blossom_node_list.size(); + + std::vector subblossoms; + _blossom_set->split(blossom, std::back_inserter(subblossoms)); + int b = _blossom_set->find(base); + int ib = -1; + for (int i = 0; i < int(subblossoms.size()); ++i) { + if (subblossoms[i] == b) { ib = i; break; } + } + + for (int i = 1; i < int(subblossoms.size()); i += 2) { + int sb = subblossoms[(ib + i) % subblossoms.size()]; + int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; + + Arc m = (*_blossom_data)[tb].next; + extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m)); + extractBlossom(tb, _graph.source(m), m); + } + extractBlossom(subblossoms[ib], base, matching); + + int en = _blossom_node_list.size(); + + _blossom_potential.push_back(BlossomVariable(bn, en, pot)); + } + } + + void extractMatching() { + std::vector blossoms; + for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { + blossoms.push_back(c); + } + + for (int i = 0; i < int(blossoms.size()); ++i) { + + Value offset = (*_blossom_data)[blossoms[i]].offset; + (*_blossom_data)[blossoms[i]].pot += 2 * offset; + for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); + n != INVALID; ++n) { + (*_node_data)[(*_node_index)[n]].pot -= offset; + } + + Arc matching = (*_blossom_data)[blossoms[i]].next; + Node base = _graph.source(matching); + extractBlossom(blossoms[i], base, matching); + } + } + + public: + + /// \brief Constructor + /// + /// Constructor. + MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight) + : _graph(graph), _weight(weight), _matching(0), + _node_potential(0), _blossom_potential(), _blossom_node_list(), + _node_num(0), _blossom_num(0), + + _blossom_index(0), _blossom_set(0), _blossom_data(0), + _node_index(0), _node_heap_index(0), _node_data(0), + _tree_set_index(0), _tree_set(0), + + _delta2_index(0), _delta2(0), + _delta3_index(0), _delta3(0), + _delta4_index(0), _delta4(0), + + _delta_sum() {} + + ~MaxWeightedPerfectMatching() { + destroyStructures(); + } + + /// \name Execution Control + /// The simplest way to execute the algorithm is to use the + /// \ref run() member function. + + ///@{ + + /// \brief Initialize the algorithm + /// + /// This function initializes the algorithm. + void init() { + createStructures(); + + for (ArcIt e(_graph); e != INVALID; ++e) { + (*_node_heap_index)[e] = BinHeap::PRE_HEAP; + } + for (EdgeIt e(_graph); e != INVALID; ++e) { + (*_delta3_index)[e] = _delta3->PRE_HEAP; + } + for (int i = 0; i < _blossom_num; ++i) { + (*_delta2_index)[i] = _delta2->PRE_HEAP; + (*_delta4_index)[i] = _delta4->PRE_HEAP; + } + + int index = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + Value max = - std::numeric_limits::max(); + for (OutArcIt e(_graph, n); e != INVALID; ++e) { + if (_graph.target(e) == n) continue; + if ((dualScale * _weight[e]) / 2 > max) { + max = (dualScale * _weight[e]) / 2; + } + } + (*_node_index)[n] = index; + (*_node_data)[index].pot = max; + int blossom = + _blossom_set->insert(n, std::numeric_limits::max()); + + _tree_set->insert(blossom); + + (*_blossom_data)[blossom].status = EVEN; + (*_blossom_data)[blossom].pred = INVALID; + (*_blossom_data)[blossom].next = INVALID; + (*_blossom_data)[blossom].pot = 0; + (*_blossom_data)[blossom].offset = 0; + ++index; + } + for (EdgeIt e(_graph); e != INVALID; ++e) { + int si = (*_node_index)[_graph.u(e)]; + int ti = (*_node_index)[_graph.v(e)]; + if (_graph.u(e) != _graph.v(e)) { + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - + dualScale * _weight[e]) / 2); + } + } + } + + /// \brief Start the algorithm + /// + /// This function starts the algorithm. + /// + /// \pre \ref init() must be called before using this function. + bool start() { + enum OpType { + D2, D3, D4 + }; + + int unmatched = _node_num; + while (unmatched > 0) { + Value d2 = !_delta2->empty() ? + _delta2->prio() : std::numeric_limits::max(); + + Value d3 = !_delta3->empty() ? + _delta3->prio() : std::numeric_limits::max(); + + Value d4 = !_delta4->empty() ? + _delta4->prio() : std::numeric_limits::max(); + + _delta_sum = d2; OpType ot = D2; + if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } + if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } + + if (_delta_sum == std::numeric_limits::max()) { + return false; + } + + switch (ot) { + case D2: + { + int blossom = _delta2->top(); + Node n = _blossom_set->classTop(blossom); + Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); + extendOnArc(e); + } + break; + case D3: + { + Edge e = _delta3->top(); + + int left_blossom = _blossom_set->find(_graph.u(e)); + int right_blossom = _blossom_set->find(_graph.v(e)); + + if (left_blossom == right_blossom) { + _delta3->pop(); + } else { + int left_tree = _tree_set->find(left_blossom); + int right_tree = _tree_set->find(right_blossom); + + if (left_tree == right_tree) { + shrinkOnEdge(e, left_tree); + } else { + augmentOnEdge(e); + unmatched -= 2; + } + } + } break; + case D4: + splitBlossom(_delta4->top()); + break; + } + } + extractMatching(); + return true; + } + + /// \brief Run the algorithm. + /// + /// This method runs the \c %MaxWeightedPerfectMatching algorithm. + /// + /// \note mwpm.run() is just a shortcut of the following code. + /// \code + /// mwpm.init(); + /// mwpm.start(); + /// \endcode + bool run() { + init(); + return start(); + } + + /// @} + + /// \name Primal Solution + /// Functions to get the primal solution, i.e. the maximum weighted + /// perfect matching.\n + /// Either \ref run() or \ref start() function should be called before + /// using them. + + /// @{ + + /// \brief Return the weight of the matching. + /// + /// This function returns the weight of the found matching. + /// + /// \pre Either run() or start() must be called before using this function. + Value matchingWeight() const { + Value sum = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_matching)[n] != INVALID) { + sum += _weight[(*_matching)[n]]; + } + } + return sum /= 2; + } + + /// \brief Return \c true if the given edge is in the matching. + /// + /// This function returns \c true if the given edge is in the found + /// matching. + /// + /// \pre Either run() or start() must be called before using this function. + bool matching(const Edge& edge) const { + return static_cast((*_matching)[_graph.u(edge)]) == edge; + } + + /// \brief Return the matching arc (or edge) incident to the given node. + /// + /// This function returns the matching arc (or edge) incident to the + /// given node in the found matching or \c INVALID if the node is + /// not covered by the matching. + /// + /// \pre Either run() or start() must be called before using this function. + Arc matching(const Node& node) const { + return (*_matching)[node]; + } + + /// \brief Return a const reference to the matching map. + /// + /// This function returns a const reference to a node map that stores + /// the matching arc (or edge) incident to each node. + const MatchingMap& matchingMap() const { + return *_matching; + } + + /// \brief Return the mate of the given node. + /// + /// This function returns the mate of the given node in the found + /// matching or \c INVALID if the node is not covered by the matching. + /// + /// \pre Either run() or start() must be called before using this function. + Node mate(const Node& node) const { + return _graph.target((*_matching)[node]); + } + + /// @} + + /// \name Dual Solution + /// Functions to get the dual solution.\n + /// Either \ref run() or \ref start() function should be called before + /// using them. + + /// @{ + + /// \brief Return the value of the dual solution. + /// + /// This function returns the value of the dual solution. + /// It should be equal to the primal value scaled by \ref dualScale + /// "dual scale". + /// + /// \pre Either run() or start() must be called before using this function. + Value dualValue() const { + Value sum = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + sum += nodeValue(n); + } + for (int i = 0; i < blossomNum(); ++i) { + sum += blossomValue(i) * (blossomSize(i) / 2); + } + return sum; + } + + /// \brief Return the dual value (potential) of the given node. + /// + /// This function returns the dual value (potential) of the given node. + /// + /// \pre Either run() or start() must be called before using this function. + Value nodeValue(const Node& n) const { + return (*_node_potential)[n]; + } + + /// \brief Return the number of the blossoms in the basis. + /// + /// This function returns the number of the blossoms in the basis. + /// + /// \pre Either run() or start() must be called before using this function. + /// \see BlossomIt + int blossomNum() const { + return _blossom_potential.size(); + } + + /// \brief Return the number of the nodes in the given blossom. + /// + /// This function returns the number of the nodes in the given blossom. + /// + /// \pre Either run() or start() must be called before using this function. + /// \see BlossomIt + int blossomSize(int k) const { + return _blossom_potential[k].end - _blossom_potential[k].begin; + } + + /// \brief Return the dual value (ptential) of the given blossom. + /// + /// This function returns the dual value (ptential) of the given blossom. + /// + /// \pre Either run() or start() must be called before using this function. + Value blossomValue(int k) const { + return _blossom_potential[k].value; + } + + /// \brief Iterator for obtaining the nodes of a blossom. + /// + /// This class provides an iterator for obtaining the nodes of the + /// given blossom. It lists a subset of the nodes. + /// Before using this iterator, you must allocate a + /// MaxWeightedPerfectMatching class and execute it. + class BlossomIt { + public: + + /// \brief Constructor. + /// + /// Constructor to get the nodes of the given variable. + /// + /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" + /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" + /// must be called before initializing this iterator. + BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) + : _algorithm(&algorithm) + { + _index = _algorithm->_blossom_potential[variable].begin; + _last = _algorithm->_blossom_potential[variable].end; + } + + /// \brief Conversion to \c Node. + /// + /// Conversion to \c Node. + operator Node() const { + return _algorithm->_blossom_node_list[_index]; + } + + /// \brief Increment operator. + /// + /// Increment operator. + BlossomIt& operator++() { + ++_index; + return *this; + } + + /// \brief Validity checking + /// + /// This function checks whether the iterator is invalid. + bool operator==(Invalid) const { return _index == _last; } + + /// \brief Validity checking + /// + /// This function checks whether the iterator is valid. + bool operator!=(Invalid) const { return _index != _last; } + + private: + const MaxWeightedPerfectMatching* _algorithm; + int _last; + int _index; + }; + + /// @} + + }; + +} //END OF NAMESPACE LEMON + +#endif //LEMON_MAX_MATCHING_H diff --git a/lemon/max_matching.h b/lemon/max_matching.h deleted file mode 100644 --- a/lemon/max_matching.h +++ /dev/null @@ -1,3107 +0,0 @@ -/* -*- mode: C++; indent-tabs-mode: nil; -*- - * - * This file is a part of LEMON, a generic C++ optimization library. - * - * Copyright (C) 2003-2009 - * 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. - * - */ - -#ifndef LEMON_MAX_MATCHING_H -#define LEMON_MAX_MATCHING_H - -#include -#include -#include -#include - -#include -#include -#include -#include - -///\ingroup matching -///\file -///\brief Maximum matching algorithms in general graphs. - -namespace lemon { - - /// \ingroup matching - /// - /// \brief Edmonds' alternating forest maximum matching algorithm. - /// - /// This class implements Edmonds' alternating forest matching - /// algorithm. The algorithm can be started from an arbitrary initial - /// matching (the default is the empty one) - /// - /// The dual solution of the problem is a map of the nodes to - /// 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 factor-critical components - /// minus the number of barrier nodes is a lower bound on the - /// unmatched nodes, and the matching is optimal if and only if this bound is - /// tight. This decomposition can be attained by calling \c - /// decomposition() after running the algorithm. - /// - /// \param GR The graph type the algorithm runs on. - template - class MaxMatching { - public: - - typedef GR Graph; - typedef typename Graph::template NodeMap - MatchingMap; - - ///\brief Indicates the Gallai-Edmonds decomposition of the graph. - /// - ///Indicates the Gallai-Edmonds decomposition of the graph. 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 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)[node] = arc; - - Node n = node; - while (n != base) { - n = _graph.target((*_matching)[n]); - Arc a = (*_ear)[n]; - n = _graph.target(a); - (*_ear)[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)[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)[_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)[node] = arc; - - Node n = node; - while (n != base) { - n = _graph.target((*_matching)[n]); - Arc a = (*_ear)[n]; - n = _graph.target(a); - (*_ear)[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)[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)[_blossom_set->find(nca)] = nca; - } - - - - void extendOnArc(const Arc& a) { - Node base = _graph.source(a); - Node odd = _graph.target(a); - - (*_ear)[odd] = _graph.oppositeArc(a); - Node even = _graph.target((*_matching)[odd]); - (*_blossom_rep)[_blossom_set->insert(even)] = even; - (*_status)[odd] = ODD; - (*_status)[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)[odd] = _graph.oppositeArc(a); - (*_status)[odd] = MATCHED; - - Arc arc = (*_matching)[even]; - (*_matching)[even] = a; - - while (arc != INVALID) { - odd = _graph.target(arc); - arc = (*_ear)[odd]; - even = _graph.target(arc); - (*_matching)[odd] = arc; - arc = (*_matching)[even]; - (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]); - } - - for (typename TreeSet::ItemIt it(*_tree_set, tree); - it != INVALID; ++it) { - if ((*_status)[it] == ODD) { - (*_status)[it] = MATCHED; - } else { - int blossom = _blossom_set->find(it); - for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom); - jt != INVALID; ++jt) { - (*_status)[jt] = MATCHED; - } - _blossom_set->eraseClass(blossom); - } - } - _tree_set->eraseClass(tree); - - } - - public: - - /// \brief Constructor - /// - /// 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 - /// \c run() member function. - /// \n - - /// If you need better control on the execution, you must call - /// \ref init(), \ref greedyInit() or \ref matchingInit() - /// functions first, then you can start the algorithm with the \ref - /// startSparse() or startDense() functions. - - ///@{ - - /// \brief Sets the actual matching to the empty matching. - /// - /// Sets the actual matching to the empty matching. - /// - void init() { - createStructures(); - for(NodeIt n(_graph); n != INVALID; ++n) { - (*_matching)[n] = INVALID; - (*_status)[n] = UNMATCHED; - } - } - - ///\brief Finds an initial matching in a greedy way - /// - ///It finds an initial matching in a greedy way. - void greedyInit() { - createStructures(); - for (NodeIt n(_graph); n != INVALID; ++n) { - (*_matching)[n] = INVALID; - (*_status)[n] = UNMATCHED; - } - 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)[n] = a; - (*_status)[n] = MATCHED; - (*_matching)[v] = _graph.oppositeArc(a); - (*_status)[v] = MATCHED; - break; - } - } - } - } - } - - - /// \brief Initialize the matching from a map containing. - /// - /// 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 \c true if the map contains a matching. - template - bool matchingInit(const MatchingMap& matching) { - createStructures(); - - for (NodeIt n(_graph); n != INVALID; ++n) { - (*_matching)[n] = INVALID; - (*_status)[n] = UNMATCHED; - } - for(EdgeIt e(_graph); e!=INVALID; ++e) { - if (matching[e]) { - - Node u = _graph.u(e); - if ((*_matching)[u] != INVALID) return false; - (*_matching)[u] = _graph.direct(e, true); - (*_status)[u] = MATCHED; - - Node v = _graph.v(e); - if ((*_matching)[v] != INVALID) return false; - (*_matching)[v] = _graph.direct(e, false); - (*_status)[v] = MATCHED; - } - } - return true; - } - - /// \brief Starts Edmonds' algorithm - /// - /// 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)[n] = EVEN; - processSparse(n); - } - } - } - - /// \brief Starts Edmonds' algorithm. - /// - /// It runs Edmonds' algorithm with a heuristic of postponing - /// shrinks, therefore resulting in 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)[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(_graph) < 2 * countNodes(_graph)) { - greedyInit(); - startSparse(); - } else { - init(); - startDense(); - } - } - - /// @} - - /// \name Primal solution - /// Functions to get the primal solution, ie. the matching. - - /// @{ - - ///\brief Returns the size of the current matching. - /// - ///Returns the size of the current matching. After \ref - ///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 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 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; - } - - /// @} - - /// \name Dual solution - /// Functions to 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. - 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) const { - return (*_status)[n] == ODD; - } - - /// @} - - }; - - /// \ingroup matching - /// - /// \brief Weighted matching in general graphs - /// - /// This class provides an efficient implementation of Edmond's - /// maximum weighted 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 - /// 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[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 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 - /** \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] */ - /// - /// The algorithm can be executed with \c run() or the \c init() and - /// then the \c start() member functions. After it the matching can - /// be asked with \c matching() or mate() functions. The dual - /// solution can be get with \c nodeValue(), \c blossomNum() and \c - /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt - /// "BlossomIt" nested class, which is able to iterate on the nodes - /// of a blossom. If the value type is integral then the dual - /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4". - template > - class MaxWeightedMatching { - public: - - ///\e - typedef GR Graph; - ///\e - typedef WM WeightMap; - ///\e - typedef typename WeightMap::Value Value; - - /// \brief Scaling factor for dual solution - /// - /// Scaling factor for dual solution, it is equal to 4 or 1 - /// according to the value type. - static const int dualScale = - std::numeric_limits::is_integer ? 4 : 1; - - typedef typename Graph::template NodeMap - MatchingMap; - - private: - - TEMPLATE_GRAPH_TYPEDEFS(Graph); - - typedef typename Graph::template NodeMap NodePotential; - typedef std::vector BlossomNodeList; - - struct BlossomVariable { - int begin, end; - Value value; - - BlossomVariable(int _begin, int _end, Value _value) - : begin(_begin), end(_end), value(_value) {} - - }; - - typedef std::vector BlossomPotential; - - const Graph& _graph; - const WeightMap& _weight; - - MatchingMap* _matching; - - NodePotential* _node_potential; - - BlossomPotential _blossom_potential; - BlossomNodeList _blossom_node_list; - - int _node_num; - int _blossom_num; - - typedef RangeMap IntIntMap; - - enum Status { - EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2 - }; - - typedef HeapUnionFind BlossomSet; - struct BlossomData { - int tree; - Status status; - Arc pred, next; - Value pot, offset; - Node base; - }; - - IntNodeMap *_blossom_index; - BlossomSet *_blossom_set; - RangeMap* _blossom_data; - - IntNodeMap *_node_index; - IntArcMap *_node_heap_index; - - struct NodeData { - - NodeData(IntArcMap& node_heap_index) - : heap(node_heap_index) {} - - int blossom; - Value pot; - BinHeap heap; - std::map heap_index; - - int tree; - }; - - RangeMap* _node_data; - - typedef ExtendFindEnum TreeSet; - - IntIntMap *_tree_set_index; - TreeSet *_tree_set; - - IntNodeMap *_delta1_index; - BinHeap *_delta1; - - IntIntMap *_delta2_index; - BinHeap *_delta2; - - IntEdgeMap *_delta3_index; - BinHeap *_delta3; - - IntIntMap *_delta4_index; - BinHeap *_delta4; - - Value _delta_sum; - - void createStructures() { - _node_num = countNodes(_graph); - _blossom_num = _node_num * 3 / 2; - - if (!_matching) { - _matching = new MatchingMap(_graph); - } - if (!_node_potential) { - _node_potential = new NodePotential(_graph); - } - if (!_blossom_set) { - _blossom_index = new IntNodeMap(_graph); - _blossom_set = new BlossomSet(*_blossom_index); - _blossom_data = new RangeMap(_blossom_num); - } - - if (!_node_index) { - _node_index = new IntNodeMap(_graph); - _node_heap_index = new IntArcMap(_graph); - _node_data = new RangeMap(_node_num, - NodeData(*_node_heap_index)); - } - - if (!_tree_set) { - _tree_set_index = new IntIntMap(_blossom_num); - _tree_set = new TreeSet(*_tree_set_index); - } - if (!_delta1) { - _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 IntEdgeMap(_graph); - _delta3 = new BinHeap(*_delta3_index); - } - if (!_delta4) { - _delta4_index = new IntIntMap(_blossom_num); - _delta4 = new BinHeap(*_delta4_index); - } - } - - void destroyStructures() { - _node_num = countNodes(_graph); - _blossom_num = _node_num * 3 / 2; - - if (_matching) { - delete _matching; - } - if (_node_potential) { - delete _node_potential; - } - if (_blossom_set) { - delete _blossom_index; - delete _blossom_set; - delete _blossom_data; - } - - if (_node_index) { - delete _node_index; - delete _node_heap_index; - delete _node_data; - } - - if (_tree_set) { - delete _tree_set_index; - delete _tree_set; - } - if (_delta1) { - delete _delta1_index; - delete _delta1; - } - if (_delta2) { - delete _delta2_index; - delete _delta2; - } - if (_delta3) { - delete _delta3_index; - delete _delta3; - } - if (_delta4) { - delete _delta4_index; - delete _delta4; - } - } - - void matchedToEven(int blossom, int tree) { - if (_delta2->state(blossom) == _delta2->IN_HEAP) { - _delta2->erase(blossom); - } - - if (!_blossom_set->trivial(blossom)) { - (*_blossom_data)[blossom].pot -= - 2 * (_delta_sum - (*_blossom_data)[blossom].offset); - } - - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); - n != INVALID; ++n) { - - _blossom_set->increase(n, std::numeric_limits::max()); - int ni = (*_node_index)[n]; - - (*_node_data)[ni].heap.clear(); - (*_node_data)[ni].heap_index.clear(); - - (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset; - - _delta1->push(n, (*_node_data)[ni].pot); - - for (InArcIt e(_graph, n); e != INVALID; ++e) { - Node v = _graph.source(e); - int vb = _blossom_set->find(v); - int vi = (*_node_index)[v]; - - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - - dualScale * _weight[e]; - - if ((*_blossom_data)[vb].status == EVEN) { - if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { - _delta3->push(e, rw / 2); - } - } else if ((*_blossom_data)[vb].status == UNMATCHED) { - if (_delta3->state(e) != _delta3->IN_HEAP) { - _delta3->push(e, rw); - } - } else { - typename std::map::iterator it = - (*_node_data)[vi].heap_index.find(tree); - - if (it != (*_node_data)[vi].heap_index.end()) { - if ((*_node_data)[vi].heap[it->second] > rw) { - (*_node_data)[vi].heap.replace(it->second, e); - (*_node_data)[vi].heap.decrease(e, rw); - it->second = e; - } - } else { - (*_node_data)[vi].heap.push(e, rw); - (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); - } - - if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { - _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); - - if ((*_blossom_data)[vb].status == MATCHED) { - if (_delta2->state(vb) != _delta2->IN_HEAP) { - _delta2->push(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset){ - _delta2->decrease(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } - } - } - } - } - } - (*_blossom_data)[blossom].offset = 0; - } - - void matchedToOdd(int blossom) { - if (_delta2->state(blossom) == _delta2->IN_HEAP) { - _delta2->erase(blossom); - } - (*_blossom_data)[blossom].offset += _delta_sum; - if (!_blossom_set->trivial(blossom)) { - _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + - (*_blossom_data)[blossom].offset); - } - } - - void evenToMatched(int blossom, int tree) { - if (!_blossom_set->trivial(blossom)) { - (*_blossom_data)[blossom].pot += 2 * _delta_sum; - } - - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); - n != INVALID; ++n) { - int ni = (*_node_index)[n]; - (*_node_data)[ni].pot -= _delta_sum; - - _delta1->erase(n); - - for (InArcIt e(_graph, n); e != INVALID; ++e) { - Node v = _graph.source(e); - int vb = _blossom_set->find(v); - int vi = (*_node_index)[v]; - - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - - dualScale * _weight[e]; - - if (vb == blossom) { - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } - } else if ((*_blossom_data)[vb].status == EVEN) { - - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } - - int vt = _tree_set->find(vb); - - if (vt != tree) { - - Arc r = _graph.oppositeArc(e); - - typename std::map::iterator it = - (*_node_data)[ni].heap_index.find(vt); - - if (it != (*_node_data)[ni].heap_index.end()) { - if ((*_node_data)[ni].heap[it->second] > rw) { - (*_node_data)[ni].heap.replace(it->second, r); - (*_node_data)[ni].heap.decrease(r, rw); - it->second = r; - } - } else { - (*_node_data)[ni].heap.push(r, rw); - (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); - } - - if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { - _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); - - if (_delta2->state(blossom) != _delta2->IN_HEAP) { - _delta2->push(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); - } else if ((*_delta2)[blossom] > - _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset){ - _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); - } - } - } - - } else if ((*_blossom_data)[vb].status == UNMATCHED) { - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } - } else { - - typename std::map::iterator it = - (*_node_data)[vi].heap_index.find(tree); - - if (it != (*_node_data)[vi].heap_index.end()) { - (*_node_data)[vi].heap.erase(it->second); - (*_node_data)[vi].heap_index.erase(it); - if ((*_node_data)[vi].heap.empty()) { - _blossom_set->increase(v, std::numeric_limits::max()); - } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { - _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); - } - - if ((*_blossom_data)[vb].status == MATCHED) { - if (_blossom_set->classPrio(vb) == - std::numeric_limits::max()) { - _delta2->erase(vb); - } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset) { - _delta2->increase(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } - } - } - } - } - } - } - - void oddToMatched(int blossom) { - (*_blossom_data)[blossom].offset -= _delta_sum; - - if (_blossom_set->classPrio(blossom) != - std::numeric_limits::max()) { - _delta2->push(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); - } - - if (!_blossom_set->trivial(blossom)) { - _delta4->erase(blossom); - } - } - - void oddToEven(int blossom, int tree) { - if (!_blossom_set->trivial(blossom)) { - _delta4->erase(blossom); - (*_blossom_data)[blossom].pot -= - 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); - } - - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); - n != INVALID; ++n) { - int ni = (*_node_index)[n]; - - _blossom_set->increase(n, std::numeric_limits::max()); - - (*_node_data)[ni].heap.clear(); - (*_node_data)[ni].heap_index.clear(); - (*_node_data)[ni].pot += - 2 * _delta_sum - (*_blossom_data)[blossom].offset; - - _delta1->push(n, (*_node_data)[ni].pot); - - for (InArcIt e(_graph, n); e != INVALID; ++e) { - Node v = _graph.source(e); - int vb = _blossom_set->find(v); - int vi = (*_node_index)[v]; - - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - - dualScale * _weight[e]; - - if ((*_blossom_data)[vb].status == EVEN) { - if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { - _delta3->push(e, rw / 2); - } - } else if ((*_blossom_data)[vb].status == UNMATCHED) { - if (_delta3->state(e) != _delta3->IN_HEAP) { - _delta3->push(e, rw); - } - } else { - - typename std::map::iterator it = - (*_node_data)[vi].heap_index.find(tree); - - if (it != (*_node_data)[vi].heap_index.end()) { - if ((*_node_data)[vi].heap[it->second] > rw) { - (*_node_data)[vi].heap.replace(it->second, e); - (*_node_data)[vi].heap.decrease(e, rw); - it->second = e; - } - } else { - (*_node_data)[vi].heap.push(e, rw); - (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); - } - - if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { - _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); - - if ((*_blossom_data)[vb].status == MATCHED) { - if (_delta2->state(vb) != _delta2->IN_HEAP) { - _delta2->push(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset) { - _delta2->decrease(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } - } - } - } - } - } - (*_blossom_data)[blossom].offset = 0; - } - - - void matchedToUnmatched(int blossom) { - if (_delta2->state(blossom) == _delta2->IN_HEAP) { - _delta2->erase(blossom); - } - - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); - n != INVALID; ++n) { - int ni = (*_node_index)[n]; - - _blossom_set->increase(n, std::numeric_limits::max()); - - (*_node_data)[ni].heap.clear(); - (*_node_data)[ni].heap_index.clear(); - - for (OutArcIt e(_graph, n); e != INVALID; ++e) { - Node v = _graph.target(e); - int vb = _blossom_set->find(v); - int vi = (*_node_index)[v]; - - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - - dualScale * _weight[e]; - - if ((*_blossom_data)[vb].status == EVEN) { - if (_delta3->state(e) != _delta3->IN_HEAP) { - _delta3->push(e, rw); - } - } - } - } - } - - void unmatchedToMatched(int blossom) { - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); - n != INVALID; ++n) { - int ni = (*_node_index)[n]; - - for (InArcIt e(_graph, n); e != INVALID; ++e) { - Node v = _graph.source(e); - int vb = _blossom_set->find(v); - int vi = (*_node_index)[v]; - - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - - dualScale * _weight[e]; - - if (vb == blossom) { - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } - } else if ((*_blossom_data)[vb].status == EVEN) { - - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } - - int vt = _tree_set->find(vb); - - Arc r = _graph.oppositeArc(e); - - typename std::map::iterator it = - (*_node_data)[ni].heap_index.find(vt); - - if (it != (*_node_data)[ni].heap_index.end()) { - if ((*_node_data)[ni].heap[it->second] > rw) { - (*_node_data)[ni].heap.replace(it->second, r); - (*_node_data)[ni].heap.decrease(r, rw); - it->second = r; - } - } else { - (*_node_data)[ni].heap.push(r, rw); - (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); - } - - if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { - _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); - - if (_delta2->state(blossom) != _delta2->IN_HEAP) { - _delta2->push(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); - } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)- - (*_blossom_data)[blossom].offset){ - _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); - } - } - - } else if ((*_blossom_data)[vb].status == UNMATCHED) { - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } - } - } - } - } - - void alternatePath(int even, int tree) { - int odd; - - evenToMatched(even, tree); - (*_blossom_data)[even].status = MATCHED; - - while ((*_blossom_data)[even].pred != INVALID) { - odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred)); - (*_blossom_data)[odd].status = MATCHED; - oddToMatched(odd); - (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; - - even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred)); - (*_blossom_data)[even].status = MATCHED; - evenToMatched(even, tree); - (*_blossom_data)[even].next = - _graph.oppositeArc((*_blossom_data)[odd].pred); - } - - } - - void destroyTree(int tree) { - for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { - if ((*_blossom_data)[b].status == EVEN) { - (*_blossom_data)[b].status = MATCHED; - evenToMatched(b, tree); - } else if ((*_blossom_data)[b].status == ODD) { - (*_blossom_data)[b].status = MATCHED; - oddToMatched(b); - } - } - _tree_set->eraseClass(tree); - } - - - void unmatchNode(const Node& node) { - int blossom = _blossom_set->find(node); - int tree = _tree_set->find(blossom); - - alternatePath(blossom, tree); - destroyTree(tree); - - (*_blossom_data)[blossom].status = UNMATCHED; - (*_blossom_data)[blossom].base = node; - matchedToUnmatched(blossom); - } - - - void augmentOnEdge(const Edge& edge) { - - int left = _blossom_set->find(_graph.u(edge)); - int right = _blossom_set->find(_graph.v(edge)); - - if ((*_blossom_data)[left].status == EVEN) { - int left_tree = _tree_set->find(left); - alternatePath(left, left_tree); - destroyTree(left_tree); - } else { - (*_blossom_data)[left].status = MATCHED; - unmatchedToMatched(left); - } - - if ((*_blossom_data)[right].status == EVEN) { - int right_tree = _tree_set->find(right); - alternatePath(right, right_tree); - destroyTree(right_tree); - } else { - (*_blossom_data)[right].status = MATCHED; - unmatchedToMatched(right); - } - - (*_blossom_data)[left].next = _graph.direct(edge, true); - (*_blossom_data)[right].next = _graph.direct(edge, false); - } - - void extendOnArc(const Arc& arc) { - int base = _blossom_set->find(_graph.target(arc)); - int tree = _tree_set->find(base); - - int odd = _blossom_set->find(_graph.source(arc)); - _tree_set->insert(odd, tree); - (*_blossom_data)[odd].status = ODD; - matchedToOdd(odd); - (*_blossom_data)[odd].pred = arc; - - int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next)); - (*_blossom_data)[even].pred = (*_blossom_data)[even].next; - _tree_set->insert(even, tree); - (*_blossom_data)[even].status = EVEN; - matchedToEven(even, tree); - } - - void shrinkOnEdge(const Edge& edge, int tree) { - int nca = -1; - std::vector left_path, right_path; - - { - std::set left_set, right_set; - int left = _blossom_set->find(_graph.u(edge)); - left_path.push_back(left); - left_set.insert(left); - - int right = _blossom_set->find(_graph.v(edge)); - right_path.push_back(right); - right_set.insert(right); - - while (true) { - - if ((*_blossom_data)[left].pred == INVALID) break; - - left = - _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); - left_path.push_back(left); - left = - _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); - left_path.push_back(left); - - left_set.insert(left); - - if (right_set.find(left) != right_set.end()) { - nca = left; - break; - } - - if ((*_blossom_data)[right].pred == INVALID) break; - - right = - _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); - right_path.push_back(right); - right = - _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); - right_path.push_back(right); - - right_set.insert(right); - - if (left_set.find(right) != left_set.end()) { - nca = right; - break; - } - - } - - if (nca == -1) { - if ((*_blossom_data)[left].pred == INVALID) { - nca = right; - while (left_set.find(nca) == left_set.end()) { - nca = - _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); - right_path.push_back(nca); - nca = - _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); - right_path.push_back(nca); - } - } else { - nca = left; - while (right_set.find(nca) == right_set.end()) { - nca = - _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); - left_path.push_back(nca); - nca = - _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); - left_path.push_back(nca); - } - } - } - } - - std::vector subblossoms; - Arc prev; - - prev = _graph.direct(edge, true); - for (int i = 0; left_path[i] != nca; i += 2) { - subblossoms.push_back(left_path[i]); - (*_blossom_data)[left_path[i]].next = prev; - _tree_set->erase(left_path[i]); - - subblossoms.push_back(left_path[i + 1]); - (*_blossom_data)[left_path[i + 1]].status = EVEN; - oddToEven(left_path[i + 1], tree); - _tree_set->erase(left_path[i + 1]); - prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred); - } - - int k = 0; - while (right_path[k] != nca) ++k; - - subblossoms.push_back(nca); - (*_blossom_data)[nca].next = prev; - - for (int i = k - 2; i >= 0; i -= 2) { - subblossoms.push_back(right_path[i + 1]); - (*_blossom_data)[right_path[i + 1]].status = EVEN; - oddToEven(right_path[i + 1], tree); - _tree_set->erase(right_path[i + 1]); - - (*_blossom_data)[right_path[i + 1]].next = - (*_blossom_data)[right_path[i + 1]].pred; - - subblossoms.push_back(right_path[i]); - _tree_set->erase(right_path[i]); - } - - int surface = - _blossom_set->join(subblossoms.begin(), subblossoms.end()); - - for (int i = 0; i < int(subblossoms.size()); ++i) { - if (!_blossom_set->trivial(subblossoms[i])) { - (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; - } - (*_blossom_data)[subblossoms[i]].status = MATCHED; - } - - (*_blossom_data)[surface].pot = -2 * _delta_sum; - (*_blossom_data)[surface].offset = 0; - (*_blossom_data)[surface].status = EVEN; - (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; - (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; - - _tree_set->insert(surface, tree); - _tree_set->erase(nca); - } - - void splitBlossom(int blossom) { - Arc next = (*_blossom_data)[blossom].next; - Arc pred = (*_blossom_data)[blossom].pred; - - int tree = _tree_set->find(blossom); - - (*_blossom_data)[blossom].status = MATCHED; - oddToMatched(blossom); - if (_delta2->state(blossom) == _delta2->IN_HEAP) { - _delta2->erase(blossom); - } - - std::vector subblossoms; - _blossom_set->split(blossom, std::back_inserter(subblossoms)); - - Value offset = (*_blossom_data)[blossom].offset; - int b = _blossom_set->find(_graph.source(pred)); - int d = _blossom_set->find(_graph.source(next)); - - int ib = -1, id = -1; - for (int i = 0; i < int(subblossoms.size()); ++i) { - if (subblossoms[i] == b) ib = i; - if (subblossoms[i] == d) id = i; - - (*_blossom_data)[subblossoms[i]].offset = offset; - if (!_blossom_set->trivial(subblossoms[i])) { - (*_blossom_data)[subblossoms[i]].pot -= 2 * offset; - } - if (_blossom_set->classPrio(subblossoms[i]) != - std::numeric_limits::max()) { - _delta2->push(subblossoms[i], - _blossom_set->classPrio(subblossoms[i]) - - (*_blossom_data)[subblossoms[i]].offset); - } - } - - if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) { - for (int i = (id + 1) % subblossoms.size(); - i != ib; i = (i + 2) % subblossoms.size()) { - int sb = subblossoms[i]; - int tb = subblossoms[(i + 1) % subblossoms.size()]; - (*_blossom_data)[sb].next = - _graph.oppositeArc((*_blossom_data)[tb].next); - } - - for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { - int sb = subblossoms[i]; - int tb = subblossoms[(i + 1) % subblossoms.size()]; - int ub = subblossoms[(i + 2) % subblossoms.size()]; - - (*_blossom_data)[sb].status = ODD; - matchedToOdd(sb); - _tree_set->insert(sb, tree); - (*_blossom_data)[sb].pred = pred; - (*_blossom_data)[sb].next = - _graph.oppositeArc((*_blossom_data)[tb].next); - - pred = (*_blossom_data)[ub].next; - - (*_blossom_data)[tb].status = EVEN; - matchedToEven(tb, tree); - _tree_set->insert(tb, tree); - (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; - } - - (*_blossom_data)[subblossoms[id]].status = ODD; - matchedToOdd(subblossoms[id]); - _tree_set->insert(subblossoms[id], tree); - (*_blossom_data)[subblossoms[id]].next = next; - (*_blossom_data)[subblossoms[id]].pred = pred; - - } else { - - for (int i = (ib + 1) % subblossoms.size(); - i != id; i = (i + 2) % subblossoms.size()) { - int sb = subblossoms[i]; - int tb = subblossoms[(i + 1) % subblossoms.size()]; - (*_blossom_data)[sb].next = - _graph.oppositeArc((*_blossom_data)[tb].next); - } - - for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { - int sb = subblossoms[i]; - int tb = subblossoms[(i + 1) % subblossoms.size()]; - int ub = subblossoms[(i + 2) % subblossoms.size()]; - - (*_blossom_data)[sb].status = ODD; - matchedToOdd(sb); - _tree_set->insert(sb, tree); - (*_blossom_data)[sb].next = next; - (*_blossom_data)[sb].pred = - _graph.oppositeArc((*_blossom_data)[tb].next); - - (*_blossom_data)[tb].status = EVEN; - matchedToEven(tb, tree); - _tree_set->insert(tb, tree); - (*_blossom_data)[tb].pred = - (*_blossom_data)[tb].next = - _graph.oppositeArc((*_blossom_data)[ub].next); - next = (*_blossom_data)[ub].next; - } - - (*_blossom_data)[subblossoms[ib]].status = ODD; - matchedToOdd(subblossoms[ib]); - _tree_set->insert(subblossoms[ib], tree); - (*_blossom_data)[subblossoms[ib]].next = next; - (*_blossom_data)[subblossoms[ib]].pred = pred; - } - _tree_set->erase(blossom); - } - - void extractBlossom(int blossom, const Node& base, const Arc& matching) { - if (_blossom_set->trivial(blossom)) { - int bi = (*_node_index)[base]; - Value pot = (*_node_data)[bi].pot; - - (*_matching)[base] = matching; - _blossom_node_list.push_back(base); - (*_node_potential)[base] = pot; - } else { - - Value pot = (*_blossom_data)[blossom].pot; - int bn = _blossom_node_list.size(); - - std::vector subblossoms; - _blossom_set->split(blossom, std::back_inserter(subblossoms)); - int b = _blossom_set->find(base); - int ib = -1; - for (int i = 0; i < int(subblossoms.size()); ++i) { - if (subblossoms[i] == b) { ib = i; break; } - } - - for (int i = 1; i < int(subblossoms.size()); i += 2) { - int sb = subblossoms[(ib + i) % subblossoms.size()]; - int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; - - Arc m = (*_blossom_data)[tb].next; - extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m)); - extractBlossom(tb, _graph.source(m), m); - } - extractBlossom(subblossoms[ib], base, matching); - - int en = _blossom_node_list.size(); - - _blossom_potential.push_back(BlossomVariable(bn, en, pot)); - } - } - - void extractMatching() { - std::vector blossoms; - for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { - blossoms.push_back(c); - } - - for (int i = 0; i < int(blossoms.size()); ++i) { - if ((*_blossom_data)[blossoms[i]].status == MATCHED) { - - Value offset = (*_blossom_data)[blossoms[i]].offset; - (*_blossom_data)[blossoms[i]].pot += 2 * offset; - for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); - n != INVALID; ++n) { - (*_node_data)[(*_node_index)[n]].pot -= offset; - } - - Arc matching = (*_blossom_data)[blossoms[i]].next; - Node base = _graph.source(matching); - extractBlossom(blossoms[i], base, matching); - } else { - Node base = (*_blossom_data)[blossoms[i]].base; - extractBlossom(blossoms[i], base, INVALID); - } - } - } - - public: - - /// \brief Constructor - /// - /// Constructor. - MaxWeightedMatching(const Graph& graph, const WeightMap& weight) - : _graph(graph), _weight(weight), _matching(0), - _node_potential(0), _blossom_potential(), _blossom_node_list(), - _node_num(0), _blossom_num(0), - - _blossom_index(0), _blossom_set(0), _blossom_data(0), - _node_index(0), _node_heap_index(0), _node_data(0), - _tree_set_index(0), _tree_set(0), - - _delta1_index(0), _delta1(0), - _delta2_index(0), _delta2(0), - _delta3_index(0), _delta3(0), - _delta4_index(0), _delta4(0), - - _delta_sum() {} - - ~MaxWeightedMatching() { - destroyStructures(); - } - - /// \name Execution control - /// The simplest way to execute the algorithm is to use the - /// \c run() member function. - - ///@{ - - /// \brief Initialize the algorithm - /// - /// Initialize the algorithm - void init() { - createStructures(); - - for (ArcIt e(_graph); e != INVALID; ++e) { - (*_node_heap_index)[e] = BinHeap::PRE_HEAP; - } - for (NodeIt n(_graph); n != INVALID; ++n) { - (*_delta1_index)[n] = _delta1->PRE_HEAP; - } - for (EdgeIt e(_graph); e != INVALID; ++e) { - (*_delta3_index)[e] = _delta3->PRE_HEAP; - } - for (int i = 0; i < _blossom_num; ++i) { - (*_delta2_index)[i] = _delta2->PRE_HEAP; - (*_delta4_index)[i] = _delta4->PRE_HEAP; - } - - int index = 0; - for (NodeIt n(_graph); n != INVALID; ++n) { - Value max = 0; - for (OutArcIt e(_graph, n); e != INVALID; ++e) { - if (_graph.target(e) == n) continue; - if ((dualScale * _weight[e]) / 2 > max) { - max = (dualScale * _weight[e]) / 2; - } - } - (*_node_index)[n] = index; - (*_node_data)[index].pot = max; - _delta1->push(n, max); - int blossom = - _blossom_set->insert(n, std::numeric_limits::max()); - - _tree_set->insert(blossom); - - (*_blossom_data)[blossom].status = EVEN; - (*_blossom_data)[blossom].pred = INVALID; - (*_blossom_data)[blossom].next = INVALID; - (*_blossom_data)[blossom].pot = 0; - (*_blossom_data)[blossom].offset = 0; - ++index; - } - for (EdgeIt e(_graph); e != INVALID; ++e) { - int si = (*_node_index)[_graph.u(e)]; - int ti = (*_node_index)[_graph.v(e)]; - if (_graph.u(e) != _graph.v(e)) { - _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - - dualScale * _weight[e]) / 2); - } - } - } - - /// \brief Starts the algorithm - /// - /// Starts the algorithm - void start() { - enum OpType { - D1, D2, D3, D4 - }; - - int unmatched = _node_num; - while (unmatched > 0) { - Value d1 = !_delta1->empty() ? - _delta1->prio() : std::numeric_limits::max(); - - Value d2 = !_delta2->empty() ? - _delta2->prio() : std::numeric_limits::max(); - - Value d3 = !_delta3->empty() ? - _delta3->prio() : std::numeric_limits::max(); - - Value d4 = !_delta4->empty() ? - _delta4->prio() : std::numeric_limits::max(); - - _delta_sum = d1; OpType ot = D1; - if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } - if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } - if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } - - - switch (ot) { - case D1: - { - Node n = _delta1->top(); - unmatchNode(n); - --unmatched; - } - break; - case D2: - { - int blossom = _delta2->top(); - Node n = _blossom_set->classTop(blossom); - Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); - extendOnArc(e); - } - break; - case D3: - { - Edge e = _delta3->top(); - - int left_blossom = _blossom_set->find(_graph.u(e)); - int right_blossom = _blossom_set->find(_graph.v(e)); - - if (left_blossom == right_blossom) { - _delta3->pop(); - } else { - int left_tree; - if ((*_blossom_data)[left_blossom].status == EVEN) { - left_tree = _tree_set->find(left_blossom); - } else { - left_tree = -1; - ++unmatched; - } - int right_tree; - if ((*_blossom_data)[right_blossom].status == EVEN) { - right_tree = _tree_set->find(right_blossom); - } else { - right_tree = -1; - ++unmatched; - } - - if (left_tree == right_tree) { - shrinkOnEdge(e, left_tree); - } else { - augmentOnEdge(e); - unmatched -= 2; - } - } - } break; - case D4: - splitBlossom(_delta4->top()); - break; - } - } - extractMatching(); - } - - /// \brief Runs %MaxWeightedMatching algorithm. - /// - /// This method runs the %MaxWeightedMatching algorithm. - /// - /// \note mwm.run() is just a shortcut of the following code. - /// \code - /// mwm.init(); - /// mwm.start(); - /// \endcode - void run() { - init(); - start(); - } - - /// @} - - /// \name Primal solution - /// Functions to get the primal solution, ie. the matching. - - /// @{ - - /// \brief Returns the weight of the matching. - /// - /// Returns the weight of the matching. - Value matchingValue() const { - Value sum = 0; - for (NodeIt n(_graph); n != INVALID; ++n) { - if ((*_matching)[n] != INVALID) { - sum += _weight[(*_matching)[n]]; - } - } - return sum /= 2; - } - - /// \brief Returns the cardinality of the matching. - /// - /// 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. - /// - /// Returns the incident matching arc from given node. If the - /// node is not matched then it gives back \c INVALID. - Arc matching(const Node& node) const { - return (*_matching)[node]; - } - - /// \brief Returns the mate of the node. - /// - /// Returns the adjancent node in a mathcing arc. If the node is - /// not matched then it gives back \c INVALID. - Node mate(const Node& node) const { - return (*_matching)[node] != INVALID ? - _graph.target((*_matching)[node]) : INVALID; - } - - /// @} - - /// \name Dual solution - /// Functions to get the dual solution. - - /// @{ - - /// \brief Returns the value of the dual solution. - /// - /// Returns the value of the dual solution. It should be equal to - /// the primal value scaled by \ref dualScale "dual scale". - Value dualValue() const { - Value sum = 0; - for (NodeIt n(_graph); n != INVALID; ++n) { - sum += nodeValue(n); - } - for (int i = 0; i < blossomNum(); ++i) { - sum += blossomValue(i) * (blossomSize(i) / 2); - } - return sum; - } - - /// \brief Returns the value of the node. - /// - /// Returns the the value of the node. - Value nodeValue(const Node& n) const { - return (*_node_potential)[n]; - } - - /// \brief Returns the number of the blossoms in the basis. - /// - /// Returns the number of the blossoms in the basis. - /// \see BlossomIt - int blossomNum() const { - return _blossom_potential.size(); - } - - - /// \brief Returns the number of the nodes in the blossom. - /// - /// Returns the number of the nodes in the blossom. - int blossomSize(int k) const { - return _blossom_potential[k].end - _blossom_potential[k].begin; - } - - /// \brief Returns the value of the blossom. - /// - /// Returns the the value of the blossom. - /// \see BlossomIt - Value blossomValue(int k) const { - return _blossom_potential[k].value; - } - - /// \brief Iterator for obtaining the nodes of the blossom. - /// - /// Iterator for obtaining the nodes of the blossom. This class - /// provides a common lemon style iterator for listing a - /// subset of the nodes. - class BlossomIt { - public: - - /// \brief Constructor. - /// - /// Constructor to get the nodes of the variable. - BlossomIt(const MaxWeightedMatching& algorithm, int variable) - : _algorithm(&algorithm) - { - _index = _algorithm->_blossom_potential[variable].begin; - _last = _algorithm->_blossom_potential[variable].end; - } - - /// \brief Conversion to node. - /// - /// Conversion to node. - operator Node() const { - return _algorithm->_blossom_node_list[_index]; - } - - /// \brief Increment operator. - /// - /// Increment operator. - BlossomIt& operator++() { - ++_index; - return *this; - } - - /// \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; - int _last; - int _index; - }; - - /// @} - - }; - - /// \ingroup matching - /// - /// \brief Weighted perfect matching in general graphs - /// - /// This class provides an efficient implementation of Edmond's - /// 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 - /// 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[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 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 - /** \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] */ - /// - /// The algorithm can be executed with \c run() or the \c init() and - /// then the \c start() member functions. After it the matching can - /// be asked with \c matching() or mate() functions. The dual - /// solution can be get with \c nodeValue(), \c blossomNum() and \c - /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt - /// "BlossomIt" nested class which is able to iterate on the nodes - /// of a blossom. If the value type is integral then the dual - /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4". - template > - class MaxWeightedPerfectMatching { - public: - - typedef GR Graph; - typedef WM WeightMap; - typedef typename WeightMap::Value Value; - - /// \brief Scaling factor for dual solution - /// - /// Scaling factor for dual solution, it is equal to 4 or 1 - /// according to the value type. - static const int dualScale = - std::numeric_limits::is_integer ? 4 : 1; - - typedef typename Graph::template NodeMap - MatchingMap; - - private: - - TEMPLATE_GRAPH_TYPEDEFS(Graph); - - typedef typename Graph::template NodeMap NodePotential; - typedef std::vector BlossomNodeList; - - struct BlossomVariable { - int begin, end; - Value value; - - BlossomVariable(int _begin, int _end, Value _value) - : begin(_begin), end(_end), value(_value) {} - - }; - - typedef std::vector BlossomPotential; - - const Graph& _graph; - const WeightMap& _weight; - - MatchingMap* _matching; - - NodePotential* _node_potential; - - BlossomPotential _blossom_potential; - BlossomNodeList _blossom_node_list; - - int _node_num; - int _blossom_num; - - typedef RangeMap IntIntMap; - - enum Status { - EVEN = -1, MATCHED = 0, ODD = 1 - }; - - typedef HeapUnionFind BlossomSet; - struct BlossomData { - int tree; - Status status; - Arc pred, next; - Value pot, offset; - }; - - IntNodeMap *_blossom_index; - BlossomSet *_blossom_set; - RangeMap* _blossom_data; - - IntNodeMap *_node_index; - IntArcMap *_node_heap_index; - - struct NodeData { - - NodeData(IntArcMap& node_heap_index) - : heap(node_heap_index) {} - - int blossom; - Value pot; - BinHeap heap; - std::map heap_index; - - int tree; - }; - - RangeMap* _node_data; - - typedef ExtendFindEnum TreeSet; - - IntIntMap *_tree_set_index; - TreeSet *_tree_set; - - IntIntMap *_delta2_index; - BinHeap *_delta2; - - IntEdgeMap *_delta3_index; - BinHeap *_delta3; - - IntIntMap *_delta4_index; - BinHeap *_delta4; - - Value _delta_sum; - - void createStructures() { - _node_num = countNodes(_graph); - _blossom_num = _node_num * 3 / 2; - - if (!_matching) { - _matching = new MatchingMap(_graph); - } - if (!_node_potential) { - _node_potential = new NodePotential(_graph); - } - if (!_blossom_set) { - _blossom_index = new IntNodeMap(_graph); - _blossom_set = new BlossomSet(*_blossom_index); - _blossom_data = new RangeMap(_blossom_num); - } - - if (!_node_index) { - _node_index = new IntNodeMap(_graph); - _node_heap_index = new IntArcMap(_graph); - _node_data = new RangeMap(_node_num, - NodeData(*_node_heap_index)); - } - - if (!_tree_set) { - _tree_set_index = new IntIntMap(_blossom_num); - _tree_set = new TreeSet(*_tree_set_index); - } - if (!_delta2) { - _delta2_index = new IntIntMap(_blossom_num); - _delta2 = new BinHeap(*_delta2_index); - } - if (!_delta3) { - _delta3_index = new IntEdgeMap(_graph); - _delta3 = new BinHeap(*_delta3_index); - } - if (!_delta4) { - _delta4_index = new IntIntMap(_blossom_num); - _delta4 = new BinHeap(*_delta4_index); - } - } - - void destroyStructures() { - _node_num = countNodes(_graph); - _blossom_num = _node_num * 3 / 2; - - if (_matching) { - delete _matching; - } - if (_node_potential) { - delete _node_potential; - } - if (_blossom_set) { - delete _blossom_index; - delete _blossom_set; - delete _blossom_data; - } - - if (_node_index) { - delete _node_index; - delete _node_heap_index; - delete _node_data; - } - - if (_tree_set) { - delete _tree_set_index; - delete _tree_set; - } - if (_delta2) { - delete _delta2_index; - delete _delta2; - } - if (_delta3) { - delete _delta3_index; - delete _delta3; - } - if (_delta4) { - delete _delta4_index; - delete _delta4; - } - } - - void matchedToEven(int blossom, int tree) { - if (_delta2->state(blossom) == _delta2->IN_HEAP) { - _delta2->erase(blossom); - } - - if (!_blossom_set->trivial(blossom)) { - (*_blossom_data)[blossom].pot -= - 2 * (_delta_sum - (*_blossom_data)[blossom].offset); - } - - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); - n != INVALID; ++n) { - - _blossom_set->increase(n, std::numeric_limits::max()); - int ni = (*_node_index)[n]; - - (*_node_data)[ni].heap.clear(); - (*_node_data)[ni].heap_index.clear(); - - (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset; - - for (InArcIt e(_graph, n); e != INVALID; ++e) { - Node v = _graph.source(e); - int vb = _blossom_set->find(v); - int vi = (*_node_index)[v]; - - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - - dualScale * _weight[e]; - - if ((*_blossom_data)[vb].status == EVEN) { - if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { - _delta3->push(e, rw / 2); - } - } else { - typename std::map::iterator it = - (*_node_data)[vi].heap_index.find(tree); - - if (it != (*_node_data)[vi].heap_index.end()) { - if ((*_node_data)[vi].heap[it->second] > rw) { - (*_node_data)[vi].heap.replace(it->second, e); - (*_node_data)[vi].heap.decrease(e, rw); - it->second = e; - } - } else { - (*_node_data)[vi].heap.push(e, rw); - (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); - } - - if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { - _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); - - if ((*_blossom_data)[vb].status == MATCHED) { - if (_delta2->state(vb) != _delta2->IN_HEAP) { - _delta2->push(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset){ - _delta2->decrease(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } - } - } - } - } - } - (*_blossom_data)[blossom].offset = 0; - } - - void matchedToOdd(int blossom) { - if (_delta2->state(blossom) == _delta2->IN_HEAP) { - _delta2->erase(blossom); - } - (*_blossom_data)[blossom].offset += _delta_sum; - if (!_blossom_set->trivial(blossom)) { - _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + - (*_blossom_data)[blossom].offset); - } - } - - void evenToMatched(int blossom, int tree) { - if (!_blossom_set->trivial(blossom)) { - (*_blossom_data)[blossom].pot += 2 * _delta_sum; - } - - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); - n != INVALID; ++n) { - int ni = (*_node_index)[n]; - (*_node_data)[ni].pot -= _delta_sum; - - for (InArcIt e(_graph, n); e != INVALID; ++e) { - Node v = _graph.source(e); - int vb = _blossom_set->find(v); - int vi = (*_node_index)[v]; - - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - - dualScale * _weight[e]; - - if (vb == blossom) { - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } - } else if ((*_blossom_data)[vb].status == EVEN) { - - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } - - int vt = _tree_set->find(vb); - - if (vt != tree) { - - Arc r = _graph.oppositeArc(e); - - typename std::map::iterator it = - (*_node_data)[ni].heap_index.find(vt); - - if (it != (*_node_data)[ni].heap_index.end()) { - if ((*_node_data)[ni].heap[it->second] > rw) { - (*_node_data)[ni].heap.replace(it->second, r); - (*_node_data)[ni].heap.decrease(r, rw); - it->second = r; - } - } else { - (*_node_data)[ni].heap.push(r, rw); - (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); - } - - if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { - _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); - - if (_delta2->state(blossom) != _delta2->IN_HEAP) { - _delta2->push(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); - } else if ((*_delta2)[blossom] > - _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset){ - _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); - } - } - } - } else { - - typename std::map::iterator it = - (*_node_data)[vi].heap_index.find(tree); - - if (it != (*_node_data)[vi].heap_index.end()) { - (*_node_data)[vi].heap.erase(it->second); - (*_node_data)[vi].heap_index.erase(it); - if ((*_node_data)[vi].heap.empty()) { - _blossom_set->increase(v, std::numeric_limits::max()); - } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { - _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); - } - - if ((*_blossom_data)[vb].status == MATCHED) { - if (_blossom_set->classPrio(vb) == - std::numeric_limits::max()) { - _delta2->erase(vb); - } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset) { - _delta2->increase(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } - } - } - } - } - } - } - - void oddToMatched(int blossom) { - (*_blossom_data)[blossom].offset -= _delta_sum; - - if (_blossom_set->classPrio(blossom) != - std::numeric_limits::max()) { - _delta2->push(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); - } - - if (!_blossom_set->trivial(blossom)) { - _delta4->erase(blossom); - } - } - - void oddToEven(int blossom, int tree) { - if (!_blossom_set->trivial(blossom)) { - _delta4->erase(blossom); - (*_blossom_data)[blossom].pot -= - 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); - } - - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); - n != INVALID; ++n) { - int ni = (*_node_index)[n]; - - _blossom_set->increase(n, std::numeric_limits::max()); - - (*_node_data)[ni].heap.clear(); - (*_node_data)[ni].heap_index.clear(); - (*_node_data)[ni].pot += - 2 * _delta_sum - (*_blossom_data)[blossom].offset; - - for (InArcIt e(_graph, n); e != INVALID; ++e) { - Node v = _graph.source(e); - int vb = _blossom_set->find(v); - int vi = (*_node_index)[v]; - - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - - dualScale * _weight[e]; - - if ((*_blossom_data)[vb].status == EVEN) { - if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { - _delta3->push(e, rw / 2); - } - } else { - - typename std::map::iterator it = - (*_node_data)[vi].heap_index.find(tree); - - if (it != (*_node_data)[vi].heap_index.end()) { - if ((*_node_data)[vi].heap[it->second] > rw) { - (*_node_data)[vi].heap.replace(it->second, e); - (*_node_data)[vi].heap.decrease(e, rw); - it->second = e; - } - } else { - (*_node_data)[vi].heap.push(e, rw); - (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); - } - - if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { - _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); - - if ((*_blossom_data)[vb].status == MATCHED) { - if (_delta2->state(vb) != _delta2->IN_HEAP) { - _delta2->push(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset) { - _delta2->decrease(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } - } - } - } - } - } - (*_blossom_data)[blossom].offset = 0; - } - - void alternatePath(int even, int tree) { - int odd; - - evenToMatched(even, tree); - (*_blossom_data)[even].status = MATCHED; - - while ((*_blossom_data)[even].pred != INVALID) { - odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred)); - (*_blossom_data)[odd].status = MATCHED; - oddToMatched(odd); - (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; - - even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred)); - (*_blossom_data)[even].status = MATCHED; - evenToMatched(even, tree); - (*_blossom_data)[even].next = - _graph.oppositeArc((*_blossom_data)[odd].pred); - } - - } - - void destroyTree(int tree) { - for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { - if ((*_blossom_data)[b].status == EVEN) { - (*_blossom_data)[b].status = MATCHED; - evenToMatched(b, tree); - } else if ((*_blossom_data)[b].status == ODD) { - (*_blossom_data)[b].status = MATCHED; - oddToMatched(b); - } - } - _tree_set->eraseClass(tree); - } - - 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); - destroyTree(left_tree); - - int right_tree = _tree_set->find(right); - alternatePath(right, right_tree); - destroyTree(right_tree); - - (*_blossom_data)[left].next = _graph.direct(edge, true); - (*_blossom_data)[right].next = _graph.direct(edge, false); - } - - void extendOnArc(const Arc& arc) { - int base = _blossom_set->find(_graph.target(arc)); - int tree = _tree_set->find(base); - - int odd = _blossom_set->find(_graph.source(arc)); - _tree_set->insert(odd, tree); - (*_blossom_data)[odd].status = ODD; - matchedToOdd(odd); - (*_blossom_data)[odd].pred = arc; - - int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next)); - (*_blossom_data)[even].pred = (*_blossom_data)[even].next; - _tree_set->insert(even, tree); - (*_blossom_data)[even].status = EVEN; - matchedToEven(even, tree); - } - - void shrinkOnEdge(const Edge& edge, int tree) { - int nca = -1; - std::vector left_path, right_path; - - { - std::set left_set, right_set; - int left = _blossom_set->find(_graph.u(edge)); - left_path.push_back(left); - left_set.insert(left); - - int right = _blossom_set->find(_graph.v(edge)); - right_path.push_back(right); - right_set.insert(right); - - while (true) { - - if ((*_blossom_data)[left].pred == INVALID) break; - - left = - _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); - left_path.push_back(left); - left = - _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); - left_path.push_back(left); - - left_set.insert(left); - - if (right_set.find(left) != right_set.end()) { - nca = left; - break; - } - - if ((*_blossom_data)[right].pred == INVALID) break; - - right = - _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); - right_path.push_back(right); - right = - _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); - right_path.push_back(right); - - right_set.insert(right); - - if (left_set.find(right) != left_set.end()) { - nca = right; - break; - } - - } - - if (nca == -1) { - if ((*_blossom_data)[left].pred == INVALID) { - nca = right; - while (left_set.find(nca) == left_set.end()) { - nca = - _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); - right_path.push_back(nca); - nca = - _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); - right_path.push_back(nca); - } - } else { - nca = left; - while (right_set.find(nca) == right_set.end()) { - nca = - _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); - left_path.push_back(nca); - nca = - _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); - left_path.push_back(nca); - } - } - } - } - - std::vector subblossoms; - Arc prev; - - prev = _graph.direct(edge, true); - for (int i = 0; left_path[i] != nca; i += 2) { - subblossoms.push_back(left_path[i]); - (*_blossom_data)[left_path[i]].next = prev; - _tree_set->erase(left_path[i]); - - subblossoms.push_back(left_path[i + 1]); - (*_blossom_data)[left_path[i + 1]].status = EVEN; - oddToEven(left_path[i + 1], tree); - _tree_set->erase(left_path[i + 1]); - prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred); - } - - int k = 0; - while (right_path[k] != nca) ++k; - - subblossoms.push_back(nca); - (*_blossom_data)[nca].next = prev; - - for (int i = k - 2; i >= 0; i -= 2) { - subblossoms.push_back(right_path[i + 1]); - (*_blossom_data)[right_path[i + 1]].status = EVEN; - oddToEven(right_path[i + 1], tree); - _tree_set->erase(right_path[i + 1]); - - (*_blossom_data)[right_path[i + 1]].next = - (*_blossom_data)[right_path[i + 1]].pred; - - subblossoms.push_back(right_path[i]); - _tree_set->erase(right_path[i]); - } - - int surface = - _blossom_set->join(subblossoms.begin(), subblossoms.end()); - - for (int i = 0; i < int(subblossoms.size()); ++i) { - if (!_blossom_set->trivial(subblossoms[i])) { - (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; - } - (*_blossom_data)[subblossoms[i]].status = MATCHED; - } - - (*_blossom_data)[surface].pot = -2 * _delta_sum; - (*_blossom_data)[surface].offset = 0; - (*_blossom_data)[surface].status = EVEN; - (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; - (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; - - _tree_set->insert(surface, tree); - _tree_set->erase(nca); - } - - void splitBlossom(int blossom) { - Arc next = (*_blossom_data)[blossom].next; - Arc pred = (*_blossom_data)[blossom].pred; - - int tree = _tree_set->find(blossom); - - (*_blossom_data)[blossom].status = MATCHED; - oddToMatched(blossom); - if (_delta2->state(blossom) == _delta2->IN_HEAP) { - _delta2->erase(blossom); - } - - std::vector subblossoms; - _blossom_set->split(blossom, std::back_inserter(subblossoms)); - - Value offset = (*_blossom_data)[blossom].offset; - int b = _blossom_set->find(_graph.source(pred)); - int d = _blossom_set->find(_graph.source(next)); - - int ib = -1, id = -1; - for (int i = 0; i < int(subblossoms.size()); ++i) { - if (subblossoms[i] == b) ib = i; - if (subblossoms[i] == d) id = i; - - (*_blossom_data)[subblossoms[i]].offset = offset; - if (!_blossom_set->trivial(subblossoms[i])) { - (*_blossom_data)[subblossoms[i]].pot -= 2 * offset; - } - if (_blossom_set->classPrio(subblossoms[i]) != - std::numeric_limits::max()) { - _delta2->push(subblossoms[i], - _blossom_set->classPrio(subblossoms[i]) - - (*_blossom_data)[subblossoms[i]].offset); - } - } - - if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) { - for (int i = (id + 1) % subblossoms.size(); - i != ib; i = (i + 2) % subblossoms.size()) { - int sb = subblossoms[i]; - int tb = subblossoms[(i + 1) % subblossoms.size()]; - (*_blossom_data)[sb].next = - _graph.oppositeArc((*_blossom_data)[tb].next); - } - - for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { - int sb = subblossoms[i]; - int tb = subblossoms[(i + 1) % subblossoms.size()]; - int ub = subblossoms[(i + 2) % subblossoms.size()]; - - (*_blossom_data)[sb].status = ODD; - matchedToOdd(sb); - _tree_set->insert(sb, tree); - (*_blossom_data)[sb].pred = pred; - (*_blossom_data)[sb].next = - _graph.oppositeArc((*_blossom_data)[tb].next); - - pred = (*_blossom_data)[ub].next; - - (*_blossom_data)[tb].status = EVEN; - matchedToEven(tb, tree); - _tree_set->insert(tb, tree); - (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; - } - - (*_blossom_data)[subblossoms[id]].status = ODD; - matchedToOdd(subblossoms[id]); - _tree_set->insert(subblossoms[id], tree); - (*_blossom_data)[subblossoms[id]].next = next; - (*_blossom_data)[subblossoms[id]].pred = pred; - - } else { - - for (int i = (ib + 1) % subblossoms.size(); - i != id; i = (i + 2) % subblossoms.size()) { - int sb = subblossoms[i]; - int tb = subblossoms[(i + 1) % subblossoms.size()]; - (*_blossom_data)[sb].next = - _graph.oppositeArc((*_blossom_data)[tb].next); - } - - for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { - int sb = subblossoms[i]; - int tb = subblossoms[(i + 1) % subblossoms.size()]; - int ub = subblossoms[(i + 2) % subblossoms.size()]; - - (*_blossom_data)[sb].status = ODD; - matchedToOdd(sb); - _tree_set->insert(sb, tree); - (*_blossom_data)[sb].next = next; - (*_blossom_data)[sb].pred = - _graph.oppositeArc((*_blossom_data)[tb].next); - - (*_blossom_data)[tb].status = EVEN; - matchedToEven(tb, tree); - _tree_set->insert(tb, tree); - (*_blossom_data)[tb].pred = - (*_blossom_data)[tb].next = - _graph.oppositeArc((*_blossom_data)[ub].next); - next = (*_blossom_data)[ub].next; - } - - (*_blossom_data)[subblossoms[ib]].status = ODD; - matchedToOdd(subblossoms[ib]); - _tree_set->insert(subblossoms[ib], tree); - (*_blossom_data)[subblossoms[ib]].next = next; - (*_blossom_data)[subblossoms[ib]].pred = pred; - } - _tree_set->erase(blossom); - } - - void extractBlossom(int blossom, const Node& base, const Arc& matching) { - if (_blossom_set->trivial(blossom)) { - int bi = (*_node_index)[base]; - Value pot = (*_node_data)[bi].pot; - - (*_matching)[base] = matching; - _blossom_node_list.push_back(base); - (*_node_potential)[base] = pot; - } else { - - Value pot = (*_blossom_data)[blossom].pot; - int bn = _blossom_node_list.size(); - - std::vector subblossoms; - _blossom_set->split(blossom, std::back_inserter(subblossoms)); - int b = _blossom_set->find(base); - int ib = -1; - for (int i = 0; i < int(subblossoms.size()); ++i) { - if (subblossoms[i] == b) { ib = i; break; } - } - - for (int i = 1; i < int(subblossoms.size()); i += 2) { - int sb = subblossoms[(ib + i) % subblossoms.size()]; - int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; - - Arc m = (*_blossom_data)[tb].next; - extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m)); - extractBlossom(tb, _graph.source(m), m); - } - extractBlossom(subblossoms[ib], base, matching); - - int en = _blossom_node_list.size(); - - _blossom_potential.push_back(BlossomVariable(bn, en, pot)); - } - } - - void extractMatching() { - std::vector blossoms; - for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { - blossoms.push_back(c); - } - - for (int i = 0; i < int(blossoms.size()); ++i) { - - Value offset = (*_blossom_data)[blossoms[i]].offset; - (*_blossom_data)[blossoms[i]].pot += 2 * offset; - for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); - n != INVALID; ++n) { - (*_node_data)[(*_node_index)[n]].pot -= offset; - } - - Arc matching = (*_blossom_data)[blossoms[i]].next; - Node base = _graph.source(matching); - extractBlossom(blossoms[i], base, matching); - } - } - - public: - - /// \brief Constructor - /// - /// Constructor. - MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight) - : _graph(graph), _weight(weight), _matching(0), - _node_potential(0), _blossom_potential(), _blossom_node_list(), - _node_num(0), _blossom_num(0), - - _blossom_index(0), _blossom_set(0), _blossom_data(0), - _node_index(0), _node_heap_index(0), _node_data(0), - _tree_set_index(0), _tree_set(0), - - _delta2_index(0), _delta2(0), - _delta3_index(0), _delta3(0), - _delta4_index(0), _delta4(0), - - _delta_sum() {} - - ~MaxWeightedPerfectMatching() { - destroyStructures(); - } - - /// \name Execution control - /// The simplest way to execute the algorithm is to use the - /// \c run() member function. - - ///@{ - - /// \brief Initialize the algorithm - /// - /// Initialize the algorithm - void init() { - createStructures(); - - for (ArcIt e(_graph); e != INVALID; ++e) { - (*_node_heap_index)[e] = BinHeap::PRE_HEAP; - } - for (EdgeIt e(_graph); e != INVALID; ++e) { - (*_delta3_index)[e] = _delta3->PRE_HEAP; - } - for (int i = 0; i < _blossom_num; ++i) { - (*_delta2_index)[i] = _delta2->PRE_HEAP; - (*_delta4_index)[i] = _delta4->PRE_HEAP; - } - - int index = 0; - for (NodeIt n(_graph); n != INVALID; ++n) { - Value max = - std::numeric_limits::max(); - for (OutArcIt e(_graph, n); e != INVALID; ++e) { - if (_graph.target(e) == n) continue; - if ((dualScale * _weight[e]) / 2 > max) { - max = (dualScale * _weight[e]) / 2; - } - } - (*_node_index)[n] = index; - (*_node_data)[index].pot = max; - int blossom = - _blossom_set->insert(n, std::numeric_limits::max()); - - _tree_set->insert(blossom); - - (*_blossom_data)[blossom].status = EVEN; - (*_blossom_data)[blossom].pred = INVALID; - (*_blossom_data)[blossom].next = INVALID; - (*_blossom_data)[blossom].pot = 0; - (*_blossom_data)[blossom].offset = 0; - ++index; - } - for (EdgeIt e(_graph); e != INVALID; ++e) { - int si = (*_node_index)[_graph.u(e)]; - int ti = (*_node_index)[_graph.v(e)]; - if (_graph.u(e) != _graph.v(e)) { - _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - - dualScale * _weight[e]) / 2); - } - } - } - - /// \brief Starts the algorithm - /// - /// Starts the algorithm - bool start() { - enum OpType { - D2, D3, D4 - }; - - int unmatched = _node_num; - while (unmatched > 0) { - Value d2 = !_delta2->empty() ? - _delta2->prio() : std::numeric_limits::max(); - - Value d3 = !_delta3->empty() ? - _delta3->prio() : std::numeric_limits::max(); - - Value d4 = !_delta4->empty() ? - _delta4->prio() : std::numeric_limits::max(); - - _delta_sum = d2; OpType ot = D2; - if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } - if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } - - if (_delta_sum == std::numeric_limits::max()) { - return false; - } - - switch (ot) { - case D2: - { - int blossom = _delta2->top(); - Node n = _blossom_set->classTop(blossom); - Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); - extendOnArc(e); - } - break; - case D3: - { - Edge e = _delta3->top(); - - int left_blossom = _blossom_set->find(_graph.u(e)); - int right_blossom = _blossom_set->find(_graph.v(e)); - - if (left_blossom == right_blossom) { - _delta3->pop(); - } else { - int left_tree = _tree_set->find(left_blossom); - int right_tree = _tree_set->find(right_blossom); - - if (left_tree == right_tree) { - shrinkOnEdge(e, left_tree); - } else { - augmentOnEdge(e); - unmatched -= 2; - } - } - } break; - case D4: - splitBlossom(_delta4->top()); - break; - } - } - extractMatching(); - return true; - } - - /// \brief Runs %MaxWeightedPerfectMatching algorithm. - /// - /// This method runs the %MaxWeightedPerfectMatching algorithm. - /// - /// \note mwm.run() is just a shortcut of the following code. - /// \code - /// mwm.init(); - /// mwm.start(); - /// \endcode - bool run() { - init(); - return start(); - } - - /// @} - - /// \name Primal solution - /// Functions to get the primal solution, ie. the matching. - - /// @{ - - /// \brief Returns the matching value. - /// - /// Returns the matching value. - Value matchingValue() const { - Value sum = 0; - for (NodeIt n(_graph); n != INVALID; ++n) { - if ((*_matching)[n] != INVALID) { - sum += _weight[(*_matching)[n]]; - } - } - return sum /= 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 static_cast((*_matching)[_graph.u(edge)]) == edge; - } - - /// \brief Returns the incident matching edge. - /// - /// Returns the incident matching arc from given edge. - Arc matching(const Node& node) const { - return (*_matching)[node]; - } - - /// \brief Returns the mate of the node. - /// - /// Returns the adjancent node in a mathcing arc. - Node mate(const Node& node) const { - return _graph.target((*_matching)[node]); - } - - /// @} - - /// \name Dual solution - /// Functions to get the dual solution. - - /// @{ - - /// \brief Returns the value of the dual solution. - /// - /// Returns the value of the dual solution. It should be equal to - /// the primal value scaled by \ref dualScale "dual scale". - Value dualValue() const { - Value sum = 0; - for (NodeIt n(_graph); n != INVALID; ++n) { - sum += nodeValue(n); - } - for (int i = 0; i < blossomNum(); ++i) { - sum += blossomValue(i) * (blossomSize(i) / 2); - } - return sum; - } - - /// \brief Returns the value of the node. - /// - /// Returns the the value of the node. - Value nodeValue(const Node& n) const { - return (*_node_potential)[n]; - } - - /// \brief Returns the number of the blossoms in the basis. - /// - /// Returns the number of the blossoms in the basis. - /// \see BlossomIt - int blossomNum() const { - return _blossom_potential.size(); - } - - - /// \brief Returns the number of the nodes in the blossom. - /// - /// Returns the number of the nodes in the blossom. - int blossomSize(int k) const { - return _blossom_potential[k].end - _blossom_potential[k].begin; - } - - /// \brief Returns the value of the blossom. - /// - /// Returns the the value of the blossom. - /// \see BlossomIt - Value blossomValue(int k) const { - return _blossom_potential[k].value; - } - - /// \brief Iterator for obtaining the nodes of the blossom. - /// - /// Iterator for obtaining the nodes of the blossom. This class - /// provides a common lemon style iterator for listing a - /// subset of the nodes. - class BlossomIt { - public: - - /// \brief Constructor. - /// - /// Constructor to get the nodes of the variable. - BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) - : _algorithm(&algorithm) - { - _index = _algorithm->_blossom_potential[variable].begin; - _last = _algorithm->_blossom_potential[variable].end; - } - - /// \brief Conversion to node. - /// - /// Conversion to node. - operator Node() const { - return _algorithm->_blossom_node_list[_index]; - } - - /// \brief Increment operator. - /// - /// Increment operator. - BlossomIt& operator++() { - ++_index; - return *this; - } - - /// \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; - int _last; - int _index; - }; - - /// @} - - }; - - -} //END OF NAMESPACE LEMON - -#endif //LEMON_MAX_MATCHING_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -29,7 +29,7 @@ heap_test kruskal_test maps_test - max_matching_test + matching_test min_cost_arborescence_test path_test preflow_test diff --git a/test/Makefile.am b/test/Makefile.am --- a/test/Makefile.am +++ b/test/Makefile.am @@ -25,7 +25,7 @@ test/heap_test \ test/kruskal_test \ test/maps_test \ - test/max_matching_test \ + test/matching_test \ test/min_cost_arborescence_test \ test/path_test \ test/preflow_test \ @@ -70,7 +70,7 @@ test_lp_test_SOURCES = test/lp_test.cc test_maps_test_SOURCES = test/maps_test.cc test_mip_test_SOURCES = test/mip_test.cc -test_max_matching_test_SOURCES = test/max_matching_test.cc +test_matching_test_SOURCES = test/matching_test.cc test_min_cost_arborescence_test_SOURCES = test/min_cost_arborescence_test.cc test_path_test_SOURCES = test/path_test.cc test_preflow_test_SOURCES = test/preflow_test.cc diff --git a/test/euler_test.cc b/test/euler_test.cc --- a/test/euler_test.cc +++ b/test/euler_test.cc @@ -18,136 +18,206 @@ #include #include -#include +#include +#include "test_tools.h" using namespace lemon; template -void checkDiEulerIt(const Digraph& g) +void checkDiEulerIt(const Digraph& g, + const typename Digraph::Node& start = INVALID) { typename Digraph::template ArcMap visitationNumber(g, 0); - DiEulerIt e(g); + DiEulerIt e(g, start); + if (e == INVALID) return; typename Digraph::Node firstNode = g.source(e); typename Digraph::Node lastNode = g.target(e); + if (start != INVALID) { + check(firstNode == start, "checkDiEulerIt: Wrong first node"); + } - for (; e != INVALID; ++e) - { - if (e != INVALID) - { - lastNode = g.target(e); - } + for (; e != INVALID; ++e) { + if (e != INVALID) lastNode = g.target(e); ++visitationNumber[e]; } check(firstNode == lastNode, - "checkDiEulerIt: first and last node are not the same"); + "checkDiEulerIt: First and last nodes are not the same"); for (typename Digraph::ArcIt a(g); a != INVALID; ++a) { check(visitationNumber[a] == 1, - "checkDiEulerIt: not visited or multiple times visited arc found"); + "checkDiEulerIt: Not visited or multiple times visited arc found"); } } template -void checkEulerIt(const Graph& g) +void checkEulerIt(const Graph& g, + const typename Graph::Node& start = INVALID) { typename Graph::template EdgeMap visitationNumber(g, 0); - EulerIt e(g); - typename Graph::Node firstNode = g.u(e); - typename Graph::Node lastNode = g.v(e); + EulerIt e(g, start); + if (e == INVALID) return; + typename Graph::Node firstNode = g.source(typename Graph::Arc(e)); + typename Graph::Node lastNode = g.target(typename Graph::Arc(e)); + if (start != INVALID) { + check(firstNode == start, "checkEulerIt: Wrong first node"); + } - for (; e != INVALID; ++e) - { - if (e != INVALID) - { - lastNode = g.v(e); - } + for (; e != INVALID; ++e) { + if (e != INVALID) lastNode = g.target(typename Graph::Arc(e)); ++visitationNumber[e]; } check(firstNode == lastNode, - "checkEulerIt: first and last node are not the same"); + "checkEulerIt: First and last nodes are not the same"); for (typename Graph::EdgeIt e(g); e != INVALID; ++e) { check(visitationNumber[e] == 1, - "checkEulerIt: not visited or multiple times visited edge found"); + "checkEulerIt: Not visited or multiple times visited edge found"); } } int main() { typedef ListDigraph Digraph; - typedef ListGraph Graph; + typedef Undirector Graph; + + { + Digraph d; + Graph g(d); + + checkDiEulerIt(d); + checkDiEulerIt(g); + checkEulerIt(g); - Digraph digraphWithEulerianCircuit; + check(eulerian(d), "This graph is Eulerian"); + check(eulerian(g), "This graph is Eulerian"); + } { - Digraph& g = digraphWithEulerianCircuit; + Digraph d; + Graph g(d); + Digraph::Node n = d.addNode(); - Digraph::Node n0 = g.addNode(); - Digraph::Node n1 = g.addNode(); - Digraph::Node n2 = g.addNode(); + checkDiEulerIt(d); + checkDiEulerIt(g); + checkEulerIt(g); - g.addArc(n0, n1); - g.addArc(n1, n0); - g.addArc(n1, n2); - g.addArc(n2, n1); + check(eulerian(d), "This graph is Eulerian"); + check(eulerian(g), "This graph is Eulerian"); } + { + Digraph d; + Graph g(d); + Digraph::Node n = d.addNode(); + d.addArc(n, n); - Digraph digraphWithoutEulerianCircuit; + checkDiEulerIt(d); + checkDiEulerIt(g); + checkEulerIt(g); + + check(eulerian(d), "This graph is Eulerian"); + check(eulerian(g), "This graph is Eulerian"); + } { - Digraph& g = digraphWithoutEulerianCircuit; + Digraph d; + Graph g(d); + Digraph::Node n1 = d.addNode(); + Digraph::Node n2 = d.addNode(); + Digraph::Node n3 = d.addNode(); + + d.addArc(n1, n2); + d.addArc(n2, n1); + d.addArc(n2, n3); + d.addArc(n3, n2); - Digraph::Node n0 = g.addNode(); - Digraph::Node n1 = g.addNode(); - Digraph::Node n2 = g.addNode(); + checkDiEulerIt(d); + checkDiEulerIt(d, n2); + checkDiEulerIt(g); + checkDiEulerIt(g, n2); + checkEulerIt(g); + checkEulerIt(g, n2); - g.addArc(n0, n1); - g.addArc(n1, n0); - g.addArc(n1, n2); + check(eulerian(d), "This graph is Eulerian"); + check(eulerian(g), "This graph is Eulerian"); } + { + Digraph d; + Graph g(d); + Digraph::Node n1 = d.addNode(); + Digraph::Node n2 = d.addNode(); + Digraph::Node n3 = d.addNode(); + Digraph::Node n4 = d.addNode(); + Digraph::Node n5 = d.addNode(); + Digraph::Node n6 = d.addNode(); + + d.addArc(n1, n2); + d.addArc(n2, n4); + d.addArc(n1, n3); + d.addArc(n3, n4); + d.addArc(n4, n1); + d.addArc(n3, n5); + d.addArc(n5, n2); + d.addArc(n4, n6); + d.addArc(n2, n6); + d.addArc(n6, n1); + d.addArc(n6, n3); - Graph graphWithEulerianCircuit; + checkDiEulerIt(d); + checkDiEulerIt(d, n1); + checkDiEulerIt(d, n5); + + checkDiEulerIt(g); + checkDiEulerIt(g, n1); + checkDiEulerIt(g, n5); + checkEulerIt(g); + checkEulerIt(g, n1); + checkEulerIt(g, n5); + + check(eulerian(d), "This graph is Eulerian"); + check(eulerian(g), "This graph is Eulerian"); + } { - Graph& g = graphWithEulerianCircuit; + Digraph d; + Graph g(d); + Digraph::Node n0 = d.addNode(); + Digraph::Node n1 = d.addNode(); + Digraph::Node n2 = d.addNode(); + Digraph::Node n3 = d.addNode(); + Digraph::Node n4 = d.addNode(); + Digraph::Node n5 = d.addNode(); + + d.addArc(n1, n2); + d.addArc(n2, n3); + d.addArc(n3, n1); - Graph::Node n0 = g.addNode(); - Graph::Node n1 = g.addNode(); - Graph::Node n2 = g.addNode(); + checkDiEulerIt(d); + checkDiEulerIt(d, n2); - g.addEdge(n0, n1); - g.addEdge(n1, n2); - g.addEdge(n2, n0); + checkDiEulerIt(g); + checkDiEulerIt(g, n2); + checkEulerIt(g); + checkEulerIt(g, n2); + + check(!eulerian(d), "This graph is not Eulerian"); + check(!eulerian(g), "This graph is not Eulerian"); } + { + Digraph d; + Graph g(d); + Digraph::Node n1 = d.addNode(); + Digraph::Node n2 = d.addNode(); + Digraph::Node n3 = d.addNode(); + + d.addArc(n1, n2); + d.addArc(n2, n3); - Graph graphWithoutEulerianCircuit; - { - Graph& g = graphWithoutEulerianCircuit; - - Graph::Node n0 = g.addNode(); - Graph::Node n1 = g.addNode(); - Graph::Node n2 = g.addNode(); - - g.addEdge(n0, n1); - g.addEdge(n1, n2); + check(!eulerian(d), "This graph is not Eulerian"); + check(!eulerian(g), "This graph is not Eulerian"); } - checkDiEulerIt(digraphWithEulerianCircuit); - - checkEulerIt(graphWithEulerianCircuit); - - check(eulerian(digraphWithEulerianCircuit), - "this graph should have an Eulerian circuit"); - check(!eulerian(digraphWithoutEulerianCircuit), - "this graph should not have an Eulerian circuit"); - - check(eulerian(graphWithEulerianCircuit), - "this graph should have an Eulerian circuit"); - check(!eulerian(graphWithoutEulerianCircuit), - "this graph should have an Eulerian circuit"); - return 0; } diff --git a/test/matching_test.cc b/test/matching_test.cc new file mode 100644 --- /dev/null +++ b/test/matching_test.cc @@ -0,0 +1,424 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2009 + * 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 +#include +#include +#include +#include + +#include "test_tools.h" + +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 checkMaxMatchingCompile() +{ + typedef concepts::Graph Graph; + typedef Graph::Node Node; + typedef Graph::Edge Edge; + typedef Graph::EdgeMap MatMap; + + Graph g; + Node n; + Edge e; + MatMap mat(g); + + MaxMatching mat_test(g); + const MaxMatching& + const_mat_test = mat_test; + + mat_test.init(); + mat_test.greedyInit(); + mat_test.matchingInit(mat); + mat_test.startSparse(); + mat_test.startDense(); + mat_test.run(); + + const_mat_test.matchingSize(); + const_mat_test.matching(e); + const_mat_test.matching(n); + const MaxMatching::MatchingMap& mmap = + const_mat_test.matchingMap(); + e = mmap[n]; + const_mat_test.mate(n); + + MaxMatching::Status stat = + const_mat_test.status(n); + const MaxMatching::StatusMap& smap = + const_mat_test.statusMap(); + stat = smap[n]; + const_mat_test.barrier(n); +} + +void checkMaxWeightedMatchingCompile() +{ + typedef concepts::Graph Graph; + typedef Graph::Node Node; + typedef Graph::Edge Edge; + typedef Graph::EdgeMap WeightMap; + + Graph g; + Node n; + Edge e; + WeightMap w(g); + + MaxWeightedMatching mat_test(g, w); + const MaxWeightedMatching& + const_mat_test = mat_test; + + mat_test.init(); + mat_test.start(); + mat_test.run(); + + const_mat_test.matchingWeight(); + const_mat_test.matchingSize(); + const_mat_test.matching(e); + const_mat_test.matching(n); + const MaxWeightedMatching::MatchingMap& mmap = + const_mat_test.matchingMap(); + e = mmap[n]; + const_mat_test.mate(n); + + int k = 0; + const_mat_test.dualValue(); + const_mat_test.nodeValue(n); + const_mat_test.blossomNum(); + const_mat_test.blossomSize(k); + const_mat_test.blossomValue(k); +} + +void checkMaxWeightedPerfectMatchingCompile() +{ + typedef concepts::Graph Graph; + typedef Graph::Node Node; + typedef Graph::Edge Edge; + typedef Graph::EdgeMap WeightMap; + + Graph g; + Node n; + Edge e; + WeightMap w(g); + + MaxWeightedPerfectMatching mat_test(g, w); + const MaxWeightedPerfectMatching& + const_mat_test = mat_test; + + mat_test.init(); + mat_test.start(); + mat_test.run(); + + const_mat_test.matchingWeight(); + const_mat_test.matching(e); + const_mat_test.matching(n); + const MaxWeightedPerfectMatching::MatchingMap& mmap = + const_mat_test.matchingMap(); + e = mmap[n]; + const_mat_test.mate(n); + + int k = 0; + const_mat_test.dualValue(); + const_mat_test.nodeValue(n); + const_mat_test.blossomNum(); + const_mat_test.blossomSize(k); + const_mat_test.blossomValue(k); +} + +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.status(n) == MaxMatching::EVEN || + mm.matching(n) != INVALID, "Wrong Gallai-Edmonds decomposition"); + if (mm.status(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.status(graph.u(e)) != MaxMatching::EVEN || + mm.status(graph.v(e)) != MaxMatching::MATCHED, + "Wrong Gallai-Edmonds decomposition"); + + check(mm.status(graph.v(e)) != MaxMatching::EVEN || + mm.status(graph.u(e)) != MaxMatching::MATCHED, + "Wrong Gallai-Edmonds decomposition"); + + if (mm.status(graph.u(e)) != MaxMatching::ODD && + mm.status(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.status(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() { + + for (int i = 0; i < lgfn; ++i) { + SmartGraph graph; + SmartGraph::EdgeMap weight(graph); + + istringstream lgfs(lgf[i]); + graphReader(graph, lgfs). + edgeMap("weight", weight).run(); + + MaxMatching mm(graph); + mm.run(); + checkMatching(graph, mm); + + MaxWeightedMatching mwm(graph, weight); + mwm.run(); + checkWeightedMatching(graph, weight, mwm); + + MaxWeightedPerfectMatching mwpm(graph, weight); + bool perfect = mwpm.run(); + + check(perfect == (mm.matchingSize() * 2 == countNodes(graph)), + "Perfect matching found"); + + if (perfect) { + checkWeightedPerfectMatching(graph, weight, mwpm); + } + } + + return 0; +} diff --git a/test/max_matching_test.cc b/test/max_matching_test.cc deleted file mode 100644 --- a/test/max_matching_test.cc +++ /dev/null @@ -1,310 +0,0 @@ -/* -*- mode: C++; indent-tabs-mode: nil; -*- - * - * This file is a part of LEMON, a generic C++ optimization library. - * - * Copyright (C) 2003-2009 - * 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 -#include -#include - -#include "test_tools.h" - -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() { - - for (int i = 0; i < lgfn; ++i) { - SmartGraph graph; - SmartGraph::EdgeMap weight(graph); - - istringstream lgfs(lgf[i]); - graphReader(graph, lgfs). - edgeMap("weight", weight).run(); - - MaxMatching mm(graph); - mm.run(); - checkMatching(graph, mm); - - MaxWeightedMatching mwm(graph, weight); - mwm.run(); - checkWeightedMatching(graph, weight, mwm); - - MaxWeightedPerfectMatching mwpm(graph, weight); - bool perfect = mwpm.run(); - - check(perfect == (mm.matchingSize() * 2 == countNodes(graph)), - "Perfect matching found"); - - if (perfect) { - checkWeightedPerfectMatching(graph, weight, mwpm); - } - } - - return 0; -} diff --git a/tools/dimacs-solver.cc b/tools/dimacs-solver.cc --- a/tools/dimacs-solver.cc +++ b/tools/dimacs-solver.cc @@ -42,7 +42,7 @@ #include #include -#include +#include using namespace lemon; typedef SmartDigraph Digraph;