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