1.1 --- a/lemon/max_matching.h Mon Oct 13 13:56:00 2008 +0200
1.2 +++ b/lemon/max_matching.h Mon Oct 13 14:00:11 2008 +0200
1.3 @@ -31,76 +31,405 @@
1.4
1.5 ///\ingroup matching
1.6 ///\file
1.7 -///\brief Maximum matching algorithms in graph.
1.8 +///\brief Maximum matching algorithms in general graphs.
1.9
1.10 namespace lemon {
1.11
1.12 - ///\ingroup matching
1.13 + /// \ingroup matching
1.14 ///
1.15 - ///\brief Edmonds' alternating forest maximum matching algorithm.
1.16 + /// \brief Edmonds' alternating forest maximum matching algorithm.
1.17 ///
1.18 - ///This class provides Edmonds' alternating forest matching
1.19 - ///algorithm. The starting matching (if any) can be passed to the
1.20 - ///algorithm using some of init functions.
1.21 + /// This class provides Edmonds' alternating forest matching
1.22 + /// algorithm. The starting matching (if any) can be passed to the
1.23 + /// algorithm using some of init functions.
1.24 ///
1.25 - ///The dual side of a matching is a map of the nodes to
1.26 - ///MaxMatching::DecompType, having values \c D, \c A and \c C
1.27 - ///showing the Gallai-Edmonds decomposition of the digraph. The nodes
1.28 - ///in \c D induce a digraph with factor-critical components, the nodes
1.29 - ///in \c A form the barrier, and the nodes in \c C induce a digraph
1.30 - ///having a perfect matching. This decomposition can be attained by
1.31 - ///calling \c decomposition() after running the algorithm.
1.32 + /// The dual side of a matching is a map of the nodes to
1.33 + /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
1.34 + /// MATCHED/C showing the Gallai-Edmonds decomposition of the
1.35 + /// graph. The nodes in \c EVEN/D induce a graph with
1.36 + /// factor-critical components, the nodes in \c ODD/A form the
1.37 + /// barrier, and the nodes in \c MATCHED/C induce a graph having a
1.38 + /// perfect matching. The number of the fractor critical components
1.39 + /// minus the number of barrier nodes is a lower bound on the
1.40 + /// unmatched nodes, and if the matching is optimal this bound is
1.41 + /// tight. This decomposition can be attained by calling \c
1.42 + /// decomposition() after running the algorithm.
1.43 ///
1.44 - ///\param Digraph The graph type the algorithm runs on.
1.45 - template <typename Graph>
1.46 + /// \param _Graph The graph type the algorithm runs on.
1.47 + template <typename _Graph>
1.48 class MaxMatching {
1.49 -
1.50 - protected:
1.51 + public:
1.52 +
1.53 + typedef _Graph Graph;
1.54 + typedef typename Graph::template NodeMap<typename Graph::Arc>
1.55 + MatchingMap;
1.56 +
1.57 + ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
1.58 + ///
1.59 + ///Indicates the Gallai-Edmonds decomposition of the graph, which
1.60 + ///shows an upper bound on the size of a maximum matching. The
1.61 + ///nodes with Status \c EVEN/D induce a graph with factor-critical
1.62 + ///components, the nodes in \c ODD/A form the canonical barrier,
1.63 + ///and the nodes in \c MATCHED/C induce a graph having a perfect
1.64 + ///matching.
1.65 + enum Status {
1.66 + EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
1.67 + };
1.68 +
1.69 + typedef typename Graph::template NodeMap<Status> StatusMap;
1.70 +
1.71 + private:
1.72
1.73 TEMPLATE_GRAPH_TYPEDEFS(Graph);
1.74
1.75 - typedef typename Graph::template NodeMap<int> UFECrossRef;
1.76 - typedef UnionFindEnum<UFECrossRef> UFE;
1.77 - typedef std::vector<Node> NV;
1.78 -
1.79 - typedef typename Graph::template NodeMap<int> EFECrossRef;
1.80 - typedef ExtendFindEnum<EFECrossRef> EFE;
1.81 + typedef UnionFindEnum<IntNodeMap> BlossomSet;
1.82 + typedef ExtendFindEnum<IntNodeMap> TreeSet;
1.83 + typedef RangeMap<Node> NodeIntMap;
1.84 + typedef MatchingMap EarMap;
1.85 + typedef std::vector<Node> NodeQueue;
1.86 +
1.87 + const Graph& _graph;
1.88 + MatchingMap* _matching;
1.89 + StatusMap* _status;
1.90 +
1.91 + EarMap* _ear;
1.92 +
1.93 + IntNodeMap* _blossom_set_index;
1.94 + BlossomSet* _blossom_set;
1.95 + NodeIntMap* _blossom_rep;
1.96 +
1.97 + IntNodeMap* _tree_set_index;
1.98 + TreeSet* _tree_set;
1.99 +
1.100 + NodeQueue _node_queue;
1.101 + int _process, _postpone, _last;
1.102 +
1.103 + int _node_num;
1.104 +
1.105 + private:
1.106 +
1.107 + void createStructures() {
1.108 + _node_num = countNodes(_graph);
1.109 + if (!_matching) {
1.110 + _matching = new MatchingMap(_graph);
1.111 + }
1.112 + if (!_status) {
1.113 + _status = new StatusMap(_graph);
1.114 + }
1.115 + if (!_ear) {
1.116 + _ear = new EarMap(_graph);
1.117 + }
1.118 + if (!_blossom_set) {
1.119 + _blossom_set_index = new IntNodeMap(_graph);
1.120 + _blossom_set = new BlossomSet(*_blossom_set_index);
1.121 + }
1.122 + if (!_blossom_rep) {
1.123 + _blossom_rep = new NodeIntMap(_node_num);
1.124 + }
1.125 + if (!_tree_set) {
1.126 + _tree_set_index = new IntNodeMap(_graph);
1.127 + _tree_set = new TreeSet(*_tree_set_index);
1.128 + }
1.129 + _node_queue.resize(_node_num);
1.130 + }
1.131 +
1.132 + void destroyStructures() {
1.133 + if (_matching) {
1.134 + delete _matching;
1.135 + }
1.136 + if (_status) {
1.137 + delete _status;
1.138 + }
1.139 + if (_ear) {
1.140 + delete _ear;
1.141 + }
1.142 + if (_blossom_set) {
1.143 + delete _blossom_set;
1.144 + delete _blossom_set_index;
1.145 + }
1.146 + if (_blossom_rep) {
1.147 + delete _blossom_rep;
1.148 + }
1.149 + if (_tree_set) {
1.150 + delete _tree_set_index;
1.151 + delete _tree_set;
1.152 + }
1.153 + }
1.154 +
1.155 + void processDense(const Node& n) {
1.156 + _process = _postpone = _last = 0;
1.157 + _node_queue[_last++] = n;
1.158 +
1.159 + while (_process != _last) {
1.160 + Node u = _node_queue[_process++];
1.161 + for (OutArcIt a(_graph, u); a != INVALID; ++a) {
1.162 + Node v = _graph.target(a);
1.163 + if ((*_status)[v] == MATCHED) {
1.164 + extendOnArc(a);
1.165 + } else if ((*_status)[v] == UNMATCHED) {
1.166 + augmentOnArc(a);
1.167 + return;
1.168 + }
1.169 + }
1.170 + }
1.171 +
1.172 + while (_postpone != _last) {
1.173 + Node u = _node_queue[_postpone++];
1.174 +
1.175 + for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
1.176 + Node v = _graph.target(a);
1.177 +
1.178 + if ((*_status)[v] == EVEN) {
1.179 + if (_blossom_set->find(u) != _blossom_set->find(v)) {
1.180 + shrinkOnEdge(a);
1.181 + }
1.182 + }
1.183 +
1.184 + while (_process != _last) {
1.185 + Node w = _node_queue[_process++];
1.186 + for (OutArcIt b(_graph, w); b != INVALID; ++b) {
1.187 + Node x = _graph.target(b);
1.188 + if ((*_status)[x] == MATCHED) {
1.189 + extendOnArc(b);
1.190 + } else if ((*_status)[x] == UNMATCHED) {
1.191 + augmentOnArc(b);
1.192 + return;
1.193 + }
1.194 + }
1.195 + }
1.196 + }
1.197 + }
1.198 + }
1.199 +
1.200 + void processSparse(const Node& n) {
1.201 + _process = _last = 0;
1.202 + _node_queue[_last++] = n;
1.203 + while (_process != _last) {
1.204 + Node u = _node_queue[_process++];
1.205 + for (OutArcIt a(_graph, u); a != INVALID; ++a) {
1.206 + Node v = _graph.target(a);
1.207 +
1.208 + if ((*_status)[v] == EVEN) {
1.209 + if (_blossom_set->find(u) != _blossom_set->find(v)) {
1.210 + shrinkOnEdge(a);
1.211 + }
1.212 + } else if ((*_status)[v] == MATCHED) {
1.213 + extendOnArc(a);
1.214 + } else if ((*_status)[v] == UNMATCHED) {
1.215 + augmentOnArc(a);
1.216 + return;
1.217 + }
1.218 + }
1.219 + }
1.220 + }
1.221 +
1.222 + void shrinkOnEdge(const Edge& e) {
1.223 + Node nca = INVALID;
1.224 +
1.225 + {
1.226 + std::set<Node> left_set, right_set;
1.227 +
1.228 + Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
1.229 + left_set.insert(left);
1.230 +
1.231 + Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
1.232 + right_set.insert(right);
1.233 +
1.234 + while (true) {
1.235 + if ((*_matching)[left] == INVALID) break;
1.236 + left = _graph.target((*_matching)[left]);
1.237 + left = (*_blossom_rep)[_blossom_set->
1.238 + find(_graph.target((*_ear)[left]))];
1.239 + if (right_set.find(left) != right_set.end()) {
1.240 + nca = left;
1.241 + break;
1.242 + }
1.243 + left_set.insert(left);
1.244 +
1.245 + if ((*_matching)[right] == INVALID) break;
1.246 + right = _graph.target((*_matching)[right]);
1.247 + right = (*_blossom_rep)[_blossom_set->
1.248 + find(_graph.target((*_ear)[right]))];
1.249 + if (left_set.find(right) != left_set.end()) {
1.250 + nca = right;
1.251 + break;
1.252 + }
1.253 + right_set.insert(right);
1.254 + }
1.255 +
1.256 + if (nca == INVALID) {
1.257 + if ((*_matching)[left] == INVALID) {
1.258 + nca = right;
1.259 + while (left_set.find(nca) == left_set.end()) {
1.260 + nca = _graph.target((*_matching)[nca]);
1.261 + nca =(*_blossom_rep)[_blossom_set->
1.262 + find(_graph.target((*_ear)[nca]))];
1.263 + }
1.264 + } else {
1.265 + nca = left;
1.266 + while (right_set.find(nca) == right_set.end()) {
1.267 + nca = _graph.target((*_matching)[nca]);
1.268 + nca = (*_blossom_rep)[_blossom_set->
1.269 + find(_graph.target((*_ear)[nca]))];
1.270 + }
1.271 + }
1.272 + }
1.273 + }
1.274 +
1.275 + {
1.276 +
1.277 + Node node = _graph.u(e);
1.278 + Arc arc = _graph.direct(e, true);
1.279 + Node base = (*_blossom_rep)[_blossom_set->find(node)];
1.280 +
1.281 + while (base != nca) {
1.282 + _ear->set(node, arc);
1.283 +
1.284 + Node n = node;
1.285 + while (n != base) {
1.286 + n = _graph.target((*_matching)[n]);
1.287 + Arc a = (*_ear)[n];
1.288 + n = _graph.target(a);
1.289 + _ear->set(n, _graph.oppositeArc(a));
1.290 + }
1.291 + node = _graph.target((*_matching)[base]);
1.292 + _tree_set->erase(base);
1.293 + _tree_set->erase(node);
1.294 + _blossom_set->insert(node, _blossom_set->find(base));
1.295 + _status->set(node, EVEN);
1.296 + _node_queue[_last++] = node;
1.297 + arc = _graph.oppositeArc((*_ear)[node]);
1.298 + node = _graph.target((*_ear)[node]);
1.299 + base = (*_blossom_rep)[_blossom_set->find(node)];
1.300 + _blossom_set->join(_graph.target(arc), base);
1.301 + }
1.302 + }
1.303 +
1.304 + _blossom_rep->set(_blossom_set->find(nca), nca);
1.305 +
1.306 + {
1.307 +
1.308 + Node node = _graph.v(e);
1.309 + Arc arc = _graph.direct(e, false);
1.310 + Node base = (*_blossom_rep)[_blossom_set->find(node)];
1.311 +
1.312 + while (base != nca) {
1.313 + _ear->set(node, arc);
1.314 +
1.315 + Node n = node;
1.316 + while (n != base) {
1.317 + n = _graph.target((*_matching)[n]);
1.318 + Arc a = (*_ear)[n];
1.319 + n = _graph.target(a);
1.320 + _ear->set(n, _graph.oppositeArc(a));
1.321 + }
1.322 + node = _graph.target((*_matching)[base]);
1.323 + _tree_set->erase(base);
1.324 + _tree_set->erase(node);
1.325 + _blossom_set->insert(node, _blossom_set->find(base));
1.326 + _status->set(node, EVEN);
1.327 + _node_queue[_last++] = node;
1.328 + arc = _graph.oppositeArc((*_ear)[node]);
1.329 + node = _graph.target((*_ear)[node]);
1.330 + base = (*_blossom_rep)[_blossom_set->find(node)];
1.331 + _blossom_set->join(_graph.target(arc), base);
1.332 + }
1.333 + }
1.334 +
1.335 + _blossom_rep->set(_blossom_set->find(nca), nca);
1.336 + }
1.337 +
1.338 +
1.339 +
1.340 + void extendOnArc(const Arc& a) {
1.341 + Node base = _graph.source(a);
1.342 + Node odd = _graph.target(a);
1.343 +
1.344 + _ear->set(odd, _graph.oppositeArc(a));
1.345 + Node even = _graph.target((*_matching)[odd]);
1.346 + _blossom_rep->set(_blossom_set->insert(even), even);
1.347 + _status->set(odd, ODD);
1.348 + _status->set(even, EVEN);
1.349 + int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
1.350 + _tree_set->insert(odd, tree);
1.351 + _tree_set->insert(even, tree);
1.352 + _node_queue[_last++] = even;
1.353 +
1.354 + }
1.355 +
1.356 + void augmentOnArc(const Arc& a) {
1.357 + Node even = _graph.source(a);
1.358 + Node odd = _graph.target(a);
1.359 +
1.360 + int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
1.361 +
1.362 + _matching->set(odd, _graph.oppositeArc(a));
1.363 + _status->set(odd, MATCHED);
1.364 +
1.365 + Arc arc = (*_matching)[even];
1.366 + _matching->set(even, a);
1.367 +
1.368 + while (arc != INVALID) {
1.369 + odd = _graph.target(arc);
1.370 + arc = (*_ear)[odd];
1.371 + even = _graph.target(arc);
1.372 + _matching->set(odd, arc);
1.373 + arc = (*_matching)[even];
1.374 + _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
1.375 + }
1.376 +
1.377 + for (typename TreeSet::ItemIt it(*_tree_set, tree);
1.378 + it != INVALID; ++it) {
1.379 + if ((*_status)[it] == ODD) {
1.380 + _status->set(it, MATCHED);
1.381 + } else {
1.382 + int blossom = _blossom_set->find(it);
1.383 + for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
1.384 + jt != INVALID; ++jt) {
1.385 + _status->set(jt, MATCHED);
1.386 + }
1.387 + _blossom_set->eraseClass(blossom);
1.388 + }
1.389 + }
1.390 + _tree_set->eraseClass(tree);
1.391 +
1.392 + }
1.393
1.394 public:
1.395
1.396 - ///\brief Indicates the Gallai-Edmonds decomposition of the digraph.
1.397 + /// \brief Constructor
1.398 ///
1.399 - ///Indicates the Gallai-Edmonds decomposition of the digraph, which
1.400 - ///shows an upper bound on the size of a maximum matching. The
1.401 - ///nodes with DecompType \c D induce a digraph with factor-critical
1.402 - ///components, the nodes in \c A form the canonical barrier, and the
1.403 - ///nodes in \c C induce a digraph having a perfect matching.
1.404 - enum DecompType {
1.405 - D=0,
1.406 - A=1,
1.407 - C=2
1.408 - };
1.409 -
1.410 - protected:
1.411 -
1.412 - static const int HEUR_density=2;
1.413 - const Graph& g;
1.414 - typename Graph::template NodeMap<Node> _mate;
1.415 - typename Graph::template NodeMap<DecompType> position;
1.416 -
1.417 - public:
1.418 -
1.419 - MaxMatching(const Graph& _g)
1.420 - : g(_g), _mate(_g), position(_g) {}
1.421 -
1.422 - ///\brief Sets the actual matching to the empty matching.
1.423 + /// Constructor.
1.424 + MaxMatching(const Graph& graph)
1.425 + : _graph(graph), _matching(0), _status(0), _ear(0),
1.426 + _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
1.427 + _tree_set_index(0), _tree_set(0) {}
1.428 +
1.429 + ~MaxMatching() {
1.430 + destroyStructures();
1.431 + }
1.432 +
1.433 + /// \name Execution control
1.434 + /// The simplest way to execute the algorithm is to use the member
1.435 + /// \c run() member function.
1.436 + /// \n
1.437 +
1.438 + /// If you need more control on the execution, first you must call
1.439 + /// \ref init(), \ref greedyInit() or \ref matchingInit()
1.440 + /// functions, then you can start the algorithm with the \ref
1.441 + /// startParse() or startDense() functions.
1.442 +
1.443 + ///@{
1.444 +
1.445 + /// \brief Sets the actual matching to the empty matching.
1.446 ///
1.447 - ///Sets the actual matching to the empty matching.
1.448 + /// Sets the actual matching to the empty matching.
1.449 ///
1.450 void init() {
1.451 - for(NodeIt v(g); v!=INVALID; ++v) {
1.452 - _mate.set(v,INVALID);
1.453 - position.set(v,C);
1.454 + createStructures();
1.455 + for(NodeIt n(_graph); n != INVALID; ++n) {
1.456 + _matching->set(n, INVALID);
1.457 + _status->set(n, UNMATCHED);
1.458 }
1.459 }
1.460
1.461 @@ -108,88 +437,96 @@
1.462 ///
1.463 ///For initial matchig it finds a maximal greedy matching.
1.464 void greedyInit() {
1.465 - for(NodeIt v(g); v!=INVALID; ++v) {
1.466 - _mate.set(v,INVALID);
1.467 - position.set(v,C);
1.468 + createStructures();
1.469 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.470 + _matching->set(n, INVALID);
1.471 + _status->set(n, UNMATCHED);
1.472 }
1.473 - for(NodeIt v(g); v!=INVALID; ++v)
1.474 - if ( _mate[v]==INVALID ) {
1.475 - for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
1.476 - Node y=g.runningNode(e);
1.477 - if ( _mate[y]==INVALID && y!=v ) {
1.478 - _mate.set(v,y);
1.479 - _mate.set(y,v);
1.480 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.481 + if ((*_matching)[n] == INVALID) {
1.482 + for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
1.483 + Node v = _graph.target(a);
1.484 + if ((*_matching)[v] == INVALID && v != n) {
1.485 + _matching->set(n, a);
1.486 + _status->set(n, MATCHED);
1.487 + _matching->set(v, _graph.oppositeArc(a));
1.488 + _status->set(v, MATCHED);
1.489 break;
1.490 }
1.491 }
1.492 }
1.493 - }
1.494 -
1.495 - ///\brief Initialize the matching from each nodes' mate.
1.496 - ///
1.497 - ///Initialize the matching from a \c Node valued \c Node map. This
1.498 - ///map must be \e symmetric, i.e. if \c map[u]==v then \c
1.499 - ///map[v]==u must hold, and \c uv will be an arc of the initial
1.500 - ///matching.
1.501 - template <typename MateMap>
1.502 - void mateMapInit(MateMap& map) {
1.503 - for(NodeIt v(g); v!=INVALID; ++v) {
1.504 - _mate.set(v,map[v]);
1.505 - position.set(v,C);
1.506 }
1.507 }
1.508
1.509 - ///\brief Initialize the matching from a node map with the
1.510 - ///incident matching arcs.
1.511 +
1.512 + /// \brief Initialize the matching from the map containing a matching.
1.513 ///
1.514 - ///Initialize the matching from an \c Edge valued \c Node map. \c
1.515 - ///map[v] must be an \c Edge incident to \c v. This map must have
1.516 - ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c
1.517 - ///g.oppositeNode(v,map[v])==u holds, and now some arc joining \c
1.518 - ///u to \c v will be an arc of the matching.
1.519 - template<typename MatchingMap>
1.520 - void matchingMapInit(MatchingMap& map) {
1.521 - for(NodeIt v(g); v!=INVALID; ++v) {
1.522 - position.set(v,C);
1.523 - Edge e=map[v];
1.524 - if ( e!=INVALID )
1.525 - _mate.set(v,g.oppositeNode(v,e));
1.526 - else
1.527 - _mate.set(v,INVALID);
1.528 + /// Initialize the matching from a \c bool valued \c Edge map. This
1.529 + /// map must have the property that there are no two incident edges
1.530 + /// with true value, ie. it contains a matching.
1.531 + /// \return %True if the map contains a matching.
1.532 + template <typename MatchingMap>
1.533 + bool matchingInit(const MatchingMap& matching) {
1.534 + createStructures();
1.535 +
1.536 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.537 + _matching->set(n, INVALID);
1.538 + _status->set(n, UNMATCHED);
1.539 }
1.540 + for(EdgeIt e(_graph); e!=INVALID; ++e) {
1.541 + if (matching[e]) {
1.542 +
1.543 + Node u = _graph.u(e);
1.544 + if ((*_matching)[u] != INVALID) return false;
1.545 + _matching->set(u, _graph.direct(e, true));
1.546 + _status->set(u, MATCHED);
1.547 +
1.548 + Node v = _graph.v(e);
1.549 + if ((*_matching)[v] != INVALID) return false;
1.550 + _matching->set(v, _graph.direct(e, false));
1.551 + _status->set(v, MATCHED);
1.552 + }
1.553 + }
1.554 + return true;
1.555 }
1.556
1.557 - ///\brief Initialize the matching from the map containing the
1.558 - ///undirected matching arcs.
1.559 + /// \brief Starts Edmonds' algorithm
1.560 ///
1.561 - ///Initialize the matching from a \c bool valued \c Edge map. This
1.562 - ///map must have the property that there are no two incident arcs
1.563 - ///\c e, \c f with \c map[e]==map[f]==true. The arcs \c e with \c
1.564 - ///map[e]==true form the matching.
1.565 - template <typename MatchingMap>
1.566 - void matchingInit(MatchingMap& map) {
1.567 - for(NodeIt v(g); v!=INVALID; ++v) {
1.568 - _mate.set(v,INVALID);
1.569 - position.set(v,C);
1.570 - }
1.571 - for(EdgeIt e(g); e!=INVALID; ++e) {
1.572 - if ( map[e] ) {
1.573 - Node u=g.u(e);
1.574 - Node v=g.v(e);
1.575 - _mate.set(u,v);
1.576 - _mate.set(v,u);
1.577 + /// If runs the original Edmonds' algorithm.
1.578 + void startSparse() {
1.579 + for(NodeIt n(_graph); n != INVALID; ++n) {
1.580 + if ((*_status)[n] == UNMATCHED) {
1.581 + (*_blossom_rep)[_blossom_set->insert(n)] = n;
1.582 + _tree_set->insert(n);
1.583 + _status->set(n, EVEN);
1.584 + processSparse(n);
1.585 }
1.586 }
1.587 }
1.588
1.589 -
1.590 - ///\brief Runs Edmonds' algorithm.
1.591 + /// \brief Starts Edmonds' algorithm.
1.592 ///
1.593 - ///Runs Edmonds' algorithm for sparse digraphs (number of arcs <
1.594 - ///2*number of nodes), and a heuristical Edmonds' algorithm with a
1.595 - ///heuristic of postponing shrinks for dense digraphs.
1.596 + /// It runs Edmonds' algorithm with a heuristic of postponing
1.597 + /// shrinks, giving a faster algorithm for dense graphs.
1.598 + void startDense() {
1.599 + for(NodeIt n(_graph); n != INVALID; ++n) {
1.600 + if ((*_status)[n] == UNMATCHED) {
1.601 + (*_blossom_rep)[_blossom_set->insert(n)] = n;
1.602 + _tree_set->insert(n);
1.603 + _status->set(n, EVEN);
1.604 + processDense(n);
1.605 + }
1.606 + }
1.607 + }
1.608 +
1.609 +
1.610 + /// \brief Runs Edmonds' algorithm
1.611 + ///
1.612 + /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
1.613 + /// or Edmonds' algorithm with a heuristic of
1.614 + /// postponing shrinks for dense graphs.
1.615 void run() {
1.616 - if (countEdges(g) < HEUR_density * countNodes(g)) {
1.617 + if (countEdges(_graph) < 2 * countNodes(_graph)) {
1.618 greedyInit();
1.619 startSparse();
1.620 } else {
1.621 @@ -198,422 +535,75 @@
1.622 }
1.623 }
1.624
1.625 -
1.626 - ///\brief Starts Edmonds' algorithm.
1.627 - ///
1.628 - ///If runs the original Edmonds' algorithm.
1.629 - void startSparse() {
1.630 -
1.631 - typename Graph::template NodeMap<Node> ear(g,INVALID);
1.632 - //undefined for the base nodes of the blossoms (i.e. for the
1.633 - //representative elements of UFE blossom) and for the nodes in C
1.634 -
1.635 - UFECrossRef blossom_base(g);
1.636 - UFE blossom(blossom_base);
1.637 - NV rep(countNodes(g));
1.638 -
1.639 - EFECrossRef tree_base(g);
1.640 - EFE tree(tree_base);
1.641 -
1.642 - //If these UFE's would be members of the class then also
1.643 - //blossom_base and tree_base should be a member.
1.644 -
1.645 - //We build only one tree and the other vertices uncovered by the
1.646 - //matching belong to C. (They can be considered as singleton
1.647 - //trees.) If this tree can be augmented or no more
1.648 - //grow/augmentation/shrink is possible then we return to this
1.649 - //"for" cycle.
1.650 - for(NodeIt v(g); v!=INVALID; ++v) {
1.651 - if (position[v]==C && _mate[v]==INVALID) {
1.652 - rep[blossom.insert(v)] = v;
1.653 - tree.insert(v);
1.654 - position.set(v,D);
1.655 - normShrink(v, ear, blossom, rep, tree);
1.656 - }
1.657 - }
1.658 - }
1.659 -
1.660 - ///\brief Starts Edmonds' algorithm.
1.661 - ///
1.662 - ///It runs Edmonds' algorithm with a heuristic of postponing
1.663 - ///shrinks, giving a faster algorithm for dense digraphs.
1.664 - void startDense() {
1.665 -
1.666 - typename Graph::template NodeMap<Node> ear(g,INVALID);
1.667 - //undefined for the base nodes of the blossoms (i.e. for the
1.668 - //representative elements of UFE blossom) and for the nodes in C
1.669 -
1.670 - UFECrossRef blossom_base(g);
1.671 - UFE blossom(blossom_base);
1.672 - NV rep(countNodes(g));
1.673 -
1.674 - EFECrossRef tree_base(g);
1.675 - EFE tree(tree_base);
1.676 -
1.677 - //If these UFE's would be members of the class then also
1.678 - //blossom_base and tree_base should be a member.
1.679 -
1.680 - //We build only one tree and the other vertices uncovered by the
1.681 - //matching belong to C. (They can be considered as singleton
1.682 - //trees.) If this tree can be augmented or no more
1.683 - //grow/augmentation/shrink is possible then we return to this
1.684 - //"for" cycle.
1.685 - for(NodeIt v(g); v!=INVALID; ++v) {
1.686 - if ( position[v]==C && _mate[v]==INVALID ) {
1.687 - rep[blossom.insert(v)] = v;
1.688 - tree.insert(v);
1.689 - position.set(v,D);
1.690 - lateShrink(v, ear, blossom, rep, tree);
1.691 - }
1.692 - }
1.693 - }
1.694 -
1.695 -
1.696 + /// @}
1.697 +
1.698 + /// \name Primal solution
1.699 + /// Functions for get the primal solution, ie. the matching.
1.700 +
1.701 + /// @{
1.702
1.703 ///\brief Returns the size of the actual matching stored.
1.704 ///
1.705 ///Returns the size of the actual matching stored. After \ref
1.706 - ///run() it returns the size of a maximum matching in the digraph.
1.707 - int size() const {
1.708 - int s=0;
1.709 - for(NodeIt v(g); v!=INVALID; ++v) {
1.710 - if ( _mate[v]!=INVALID ) {
1.711 - ++s;
1.712 + ///run() it returns the size of the maximum matching in the graph.
1.713 + int matchingSize() const {
1.714 + int size = 0;
1.715 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.716 + if ((*_matching)[n] != INVALID) {
1.717 + ++size;
1.718 }
1.719 }
1.720 - return s/2;
1.721 + return size / 2;
1.722 }
1.723
1.724 + /// \brief Returns true when the edge is in the matching.
1.725 + ///
1.726 + /// Returns true when the edge is in the matching.
1.727 + bool matching(const Edge& edge) const {
1.728 + return edge == (*_matching)[_graph.u(edge)];
1.729 + }
1.730 +
1.731 + /// \brief Returns the matching edge incident to the given node.
1.732 + ///
1.733 + /// Returns the matching edge of a \c node in the actual matching or
1.734 + /// INVALID if the \c node is not covered by the actual matching.
1.735 + Arc matching(const Node& n) const {
1.736 + return (*_matching)[n];
1.737 + }
1.738
1.739 ///\brief Returns the mate of a node in the actual matching.
1.740 ///
1.741 - ///Returns the mate of a \c node in the actual matching.
1.742 - ///Returns INVALID if the \c node is not covered by the actual matching.
1.743 - Node mate(const Node& node) const {
1.744 - return _mate[node];
1.745 + ///Returns the mate of a \c node in the actual matching or
1.746 + ///INVALID if the \c node is not covered by the actual matching.
1.747 + Node mate(const Node& n) const {
1.748 + return (*_matching)[n] != INVALID ?
1.749 + _graph.target((*_matching)[n]) : INVALID;
1.750 }
1.751
1.752 - ///\brief Returns the matching arc incident to the given node.
1.753 - ///
1.754 - ///Returns the matching arc of a \c node in the actual matching.
1.755 - ///Returns INVALID if the \c node is not covered by the actual matching.
1.756 - Edge matchingArc(const Node& node) const {
1.757 - if (_mate[node] == INVALID) return INVALID;
1.758 - Node n = node < _mate[node] ? node : _mate[node];
1.759 - for (IncEdgeIt e(g, n); e != INVALID; ++e) {
1.760 - if (g.oppositeNode(n, e) == _mate[n]) {
1.761 - return e;
1.762 - }
1.763 - }
1.764 - return INVALID;
1.765 - }
1.766 + /// @}
1.767 +
1.768 + /// \name Dual solution
1.769 + /// Functions for get the dual solution, ie. the decomposition.
1.770 +
1.771 + /// @{
1.772
1.773 /// \brief Returns the class of the node in the Edmonds-Gallai
1.774 /// decomposition.
1.775 ///
1.776 /// Returns the class of the node in the Edmonds-Gallai
1.777 /// decomposition.
1.778 - DecompType decomposition(const Node& n) {
1.779 - return position[n] == A;
1.780 + Status decomposition(const Node& n) const {
1.781 + return (*_status)[n];
1.782 }
1.783
1.784 /// \brief Returns true when the node is in the barrier.
1.785 ///
1.786 /// Returns true when the node is in the barrier.
1.787 - bool barrier(const Node& n) {
1.788 - return position[n] == A;
1.789 + bool barrier(const Node& n) const {
1.790 + return (*_status)[n] == ODD;
1.791 }
1.792
1.793 - ///\brief Gives back the matching in a \c Node of mates.
1.794 - ///
1.795 - ///Writes the stored matching to a \c Node valued \c Node map. The
1.796 - ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
1.797 - ///map[v]==u will hold, and now \c uv is an arc of the matching.
1.798 - template <typename MateMap>
1.799 - void mateMap(MateMap& map) const {
1.800 - for(NodeIt v(g); v!=INVALID; ++v) {
1.801 - map.set(v,_mate[v]);
1.802 - }
1.803 - }
1.804 -
1.805 - ///\brief Gives back the matching in an \c Edge valued \c Node
1.806 - ///map.
1.807 - ///
1.808 - ///Writes the stored matching to an \c Edge valued \c Node
1.809 - ///map. \c map[v] will be an \c Edge incident to \c v. This
1.810 - ///map will have the property that if \c g.oppositeNode(u,map[u])
1.811 - ///== v then \c map[u]==map[v] holds, and now this arc is an arc
1.812 - ///of the matching.
1.813 - template <typename MatchingMap>
1.814 - void matchingMap(MatchingMap& map) const {
1.815 - typename Graph::template NodeMap<bool> todo(g,true);
1.816 - for(NodeIt v(g); v!=INVALID; ++v) {
1.817 - if (_mate[v]!=INVALID && v < _mate[v]) {
1.818 - Node u=_mate[v];
1.819 - for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
1.820 - if ( g.runningNode(e) == u ) {
1.821 - map.set(u,e);
1.822 - map.set(v,e);
1.823 - todo.set(u,false);
1.824 - todo.set(v,false);
1.825 - break;
1.826 - }
1.827 - }
1.828 - }
1.829 - }
1.830 - }
1.831 -
1.832 -
1.833 - ///\brief Gives back the matching in a \c bool valued \c Edge
1.834 - ///map.
1.835 - ///
1.836 - ///Writes the matching stored to a \c bool valued \c Arc
1.837 - ///map. This map will have the property that there are no two
1.838 - ///incident arcs \c e, \c f with \c map[e]==map[f]==true. The
1.839 - ///arcs \c e with \c map[e]==true form the matching.
1.840 - template<typename MatchingMap>
1.841 - void matching(MatchingMap& map) const {
1.842 - for(EdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
1.843 -
1.844 - typename Graph::template NodeMap<bool> todo(g,true);
1.845 - for(NodeIt v(g); v!=INVALID; ++v) {
1.846 - if ( todo[v] && _mate[v]!=INVALID ) {
1.847 - Node u=_mate[v];
1.848 - for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
1.849 - if ( g.runningNode(e) == u ) {
1.850 - map.set(e,true);
1.851 - todo.set(u,false);
1.852 - todo.set(v,false);
1.853 - break;
1.854 - }
1.855 - }
1.856 - }
1.857 - }
1.858 - }
1.859 -
1.860 -
1.861 - ///\brief Returns the canonical decomposition of the digraph after running
1.862 - ///the algorithm.
1.863 - ///
1.864 - ///After calling any run methods of the class, it writes the
1.865 - ///Gallai-Edmonds canonical decomposition of the digraph. \c map
1.866 - ///must be a node map of \ref DecompType 's.
1.867 - template <typename DecompositionMap>
1.868 - void decomposition(DecompositionMap& map) const {
1.869 - for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]);
1.870 - }
1.871 -
1.872 - ///\brief Returns a barrier on the nodes.
1.873 - ///
1.874 - ///After calling any run methods of the class, it writes a
1.875 - ///canonical barrier on the nodes. The odd component number of the
1.876 - ///remaining digraph minus the barrier size is a lower bound for the
1.877 - ///uncovered nodes in the digraph. The \c map must be a node map of
1.878 - ///bools.
1.879 - template <typename BarrierMap>
1.880 - void barrier(BarrierMap& barrier) {
1.881 - for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A);
1.882 - }
1.883 -
1.884 - private:
1.885 -
1.886 -
1.887 - void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,
1.888 - UFE& blossom, NV& rep, EFE& tree) {
1.889 - //We have one tree which we grow, and also shrink but only if it
1.890 - //cannot be postponed. If we augment then we return to the "for"
1.891 - //cycle of runEdmonds().
1.892 -
1.893 - std::queue<Node> Q; //queue of the totally unscanned nodes
1.894 - Q.push(v);
1.895 - std::queue<Node> R;
1.896 - //queue of the nodes which must be scanned for a possible shrink
1.897 -
1.898 - while ( !Q.empty() ) {
1.899 - Node x=Q.front();
1.900 - Q.pop();
1.901 - for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
1.902 - Node y=g.runningNode(e);
1.903 - //growOrAugment grows if y is covered by the matching and
1.904 - //augments if not. In this latter case it returns 1.
1.905 - if (position[y]==C &&
1.906 - growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
1.907 - }
1.908 - R.push(x);
1.909 - }
1.910 -
1.911 - while ( !R.empty() ) {
1.912 - Node x=R.front();
1.913 - R.pop();
1.914 -
1.915 - for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
1.916 - Node y=g.runningNode(e);
1.917 -
1.918 - if ( position[y] == D && blossom.find(x) != blossom.find(y) )
1.919 - //Recall that we have only one tree.
1.920 - shrink( x, y, ear, blossom, rep, tree, Q);
1.921 -
1.922 - while ( !Q.empty() ) {
1.923 - Node z=Q.front();
1.924 - Q.pop();
1.925 - for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
1.926 - Node w=g.runningNode(f);
1.927 - //growOrAugment grows if y is covered by the matching and
1.928 - //augments if not. In this latter case it returns 1.
1.929 - if (position[w]==C &&
1.930 - growOrAugment(w, z, ear, blossom, rep, tree, Q)) return;
1.931 - }
1.932 - R.push(z);
1.933 - }
1.934 - } //for e
1.935 - } // while ( !R.empty() )
1.936 - }
1.937 -
1.938 - void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,
1.939 - UFE& blossom, NV& rep, EFE& tree) {
1.940 - //We have one tree, which we grow and shrink. If we augment then we
1.941 - //return to the "for" cycle of runEdmonds().
1.942 -
1.943 - std::queue<Node> Q; //queue of the unscanned nodes
1.944 - Q.push(v);
1.945 - while ( !Q.empty() ) {
1.946 -
1.947 - Node x=Q.front();
1.948 - Q.pop();
1.949 -
1.950 - for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
1.951 - Node y=g.runningNode(e);
1.952 -
1.953 - switch ( position[y] ) {
1.954 - case D: //x and y must be in the same tree
1.955 - if ( blossom.find(x) != blossom.find(y))
1.956 - //x and y are in the same tree
1.957 - shrink(x, y, ear, blossom, rep, tree, Q);
1.958 - break;
1.959 - case C:
1.960 - //growOrAugment grows if y is covered by the matching and
1.961 - //augments if not. In this latter case it returns 1.
1.962 - if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
1.963 - break;
1.964 - default: break;
1.965 - }
1.966 - }
1.967 - }
1.968 - }
1.969 -
1.970 - void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear,
1.971 - UFE& blossom, NV& rep, EFE& tree,std::queue<Node>& Q) {
1.972 - //x and y are the two adjacent vertices in two blossoms.
1.973 -
1.974 - typename Graph::template NodeMap<bool> path(g,false);
1.975 -
1.976 - Node b=rep[blossom.find(x)];
1.977 - path.set(b,true);
1.978 - b=_mate[b];
1.979 - while ( b!=INVALID ) {
1.980 - b=rep[blossom.find(ear[b])];
1.981 - path.set(b,true);
1.982 - b=_mate[b];
1.983 - } //we go until the root through bases of blossoms and odd vertices
1.984 -
1.985 - Node top=y;
1.986 - Node middle=rep[blossom.find(top)];
1.987 - Node bottom=x;
1.988 - while ( !path[middle] )
1.989 - shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
1.990 - //Until we arrive to a node on the path, we update blossom, tree
1.991 - //and the positions of the odd nodes.
1.992 -
1.993 - Node base=middle;
1.994 - top=x;
1.995 - middle=rep[blossom.find(top)];
1.996 - bottom=y;
1.997 - Node blossom_base=rep[blossom.find(base)];
1.998 - while ( middle!=blossom_base )
1.999 - shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
1.1000 - //Until we arrive to a node on the path, we update blossom, tree
1.1001 - //and the positions of the odd nodes.
1.1002 -
1.1003 - rep[blossom.find(base)] = base;
1.1004 - }
1.1005 -
1.1006 - void shrinkStep(Node& top, Node& middle, Node& bottom,
1.1007 - typename Graph::template NodeMap<Node>& ear,
1.1008 - UFE& blossom, NV& rep, EFE& tree, std::queue<Node>& Q) {
1.1009 - //We traverse a blossom and update everything.
1.1010 -
1.1011 - ear.set(top,bottom);
1.1012 - Node t=top;
1.1013 - while ( t!=middle ) {
1.1014 - Node u=_mate[t];
1.1015 - t=ear[u];
1.1016 - ear.set(t,u);
1.1017 - }
1.1018 - bottom=_mate[middle];
1.1019 - position.set(bottom,D);
1.1020 - Q.push(bottom);
1.1021 - top=ear[bottom];
1.1022 - Node oldmiddle=middle;
1.1023 - middle=rep[blossom.find(top)];
1.1024 - tree.erase(bottom);
1.1025 - tree.erase(oldmiddle);
1.1026 - blossom.insert(bottom);
1.1027 - blossom.join(bottom, oldmiddle);
1.1028 - blossom.join(top, oldmiddle);
1.1029 - }
1.1030 -
1.1031 -
1.1032 -
1.1033 - bool growOrAugment(Node& y, Node& x, typename Graph::template
1.1034 - NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree,
1.1035 - std::queue<Node>& Q) {
1.1036 - //x is in a blossom in the tree, y is outside. If y is covered by
1.1037 - //the matching we grow, otherwise we augment. In this case we
1.1038 - //return 1.
1.1039 -
1.1040 - if ( _mate[y]!=INVALID ) { //grow
1.1041 - ear.set(y,x);
1.1042 - Node w=_mate[y];
1.1043 - rep[blossom.insert(w)] = w;
1.1044 - position.set(y,A);
1.1045 - position.set(w,D);
1.1046 - int t = tree.find(rep[blossom.find(x)]);
1.1047 - tree.insert(y,t);
1.1048 - tree.insert(w,t);
1.1049 - Q.push(w);
1.1050 - } else { //augment
1.1051 - augment(x, ear, blossom, rep, tree);
1.1052 - _mate.set(x,y);
1.1053 - _mate.set(y,x);
1.1054 - return true;
1.1055 - }
1.1056 - return false;
1.1057 - }
1.1058 -
1.1059 - void augment(Node x, typename Graph::template NodeMap<Node>& ear,
1.1060 - UFE& blossom, NV& rep, EFE& tree) {
1.1061 - Node v=_mate[x];
1.1062 - while ( v!=INVALID ) {
1.1063 -
1.1064 - Node u=ear[v];
1.1065 - _mate.set(v,u);
1.1066 - Node tmp=v;
1.1067 - v=_mate[u];
1.1068 - _mate.set(u,tmp);
1.1069 - }
1.1070 - int y = tree.find(rep[blossom.find(x)]);
1.1071 - for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {
1.1072 - if ( position[tit] == D ) {
1.1073 - int b = blossom.find(tit);
1.1074 - for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) {
1.1075 - position.set(bit, C);
1.1076 - }
1.1077 - blossom.eraseClass(b);
1.1078 - } else position.set(tit, C);
1.1079 - }
1.1080 - tree.eraseClass(y);
1.1081 -
1.1082 - }
1.1083 + /// @}
1.1084
1.1085 };
1.1086
1.1087 @@ -627,25 +617,28 @@
1.1088 /// \f$O(nm\log(n))\f$ time complexity.
1.1089 ///
1.1090 /// The maximum weighted matching problem is to find undirected
1.1091 - /// arcs in the digraph with maximum overall weight and no two of
1.1092 - /// them shares their endpoints. The problem can be formulated with
1.1093 - /// the next linear program:
1.1094 + /// edges in the graph with maximum overall weight and no two of
1.1095 + /// them shares their ends. The problem can be formulated with the
1.1096 + /// following linear program.
1.1097 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
1.1098 - ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
1.1099 + /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
1.1100 + \quad \forall B\in\mathcal{O}\f] */
1.1101 /// \f[x_e \ge 0\quad \forall e\in E\f]
1.1102 /// \f[\max \sum_{e\in E}x_ew_e\f]
1.1103 - /// where \f$\delta(X)\f$ is the set of arcs incident to a node in
1.1104 - /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in
1.1105 - /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
1.1106 - /// the nodes.
1.1107 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1.1108 + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
1.1109 + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
1.1110 + /// subsets of the nodes.
1.1111 ///
1.1112 /// The algorithm calculates an optimal matching and a proof of the
1.1113 /// optimality. The solution of the dual problem can be used to check
1.1114 - /// the result of the algorithm. The dual linear problem is the next:
1.1115 - /// \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]
1.1116 + /// the result of the algorithm. The dual linear problem is the
1.1117 + /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
1.1118 + z_B \ge w_{uv} \quad \forall uv\in E\f] */
1.1119 /// \f[y_u \ge 0 \quad \forall u \in V\f]
1.1120 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1.1121 - /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
1.1122 + /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1.1123 + \frac{\vert B \vert - 1}{2}z_B\f] */
1.1124 ///
1.1125 /// The algorithm can be executed with \c run() or the \c init() and
1.1126 /// then the \c start() member functions. After it the matching can
1.1127 @@ -705,16 +698,13 @@
1.1128 int _node_num;
1.1129 int _blossom_num;
1.1130
1.1131 - typedef typename Graph::template NodeMap<int> NodeIntMap;
1.1132 - typedef typename Graph::template ArcMap<int> ArcIntMap;
1.1133 - typedef typename Graph::template EdgeMap<int> EdgeIntMap;
1.1134 typedef RangeMap<int> IntIntMap;
1.1135
1.1136 enum Status {
1.1137 EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
1.1138 };
1.1139
1.1140 - typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
1.1141 + typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
1.1142 struct BlossomData {
1.1143 int tree;
1.1144 Status status;
1.1145 @@ -723,21 +713,21 @@
1.1146 Node base;
1.1147 };
1.1148
1.1149 - NodeIntMap *_blossom_index;
1.1150 + IntNodeMap *_blossom_index;
1.1151 BlossomSet *_blossom_set;
1.1152 RangeMap<BlossomData>* _blossom_data;
1.1153
1.1154 - NodeIntMap *_node_index;
1.1155 - ArcIntMap *_node_heap_index;
1.1156 + IntNodeMap *_node_index;
1.1157 + IntArcMap *_node_heap_index;
1.1158
1.1159 struct NodeData {
1.1160
1.1161 - NodeData(ArcIntMap& node_heap_index)
1.1162 + NodeData(IntArcMap& node_heap_index)
1.1163 : heap(node_heap_index) {}
1.1164
1.1165 int blossom;
1.1166 Value pot;
1.1167 - BinHeap<Value, ArcIntMap> heap;
1.1168 + BinHeap<Value, IntArcMap> heap;
1.1169 std::map<int, Arc> heap_index;
1.1170
1.1171 int tree;
1.1172 @@ -750,14 +740,14 @@
1.1173 IntIntMap *_tree_set_index;
1.1174 TreeSet *_tree_set;
1.1175
1.1176 - NodeIntMap *_delta1_index;
1.1177 - BinHeap<Value, NodeIntMap> *_delta1;
1.1178 + IntNodeMap *_delta1_index;
1.1179 + BinHeap<Value, IntNodeMap> *_delta1;
1.1180
1.1181 IntIntMap *_delta2_index;
1.1182 BinHeap<Value, IntIntMap> *_delta2;
1.1183
1.1184 - EdgeIntMap *_delta3_index;
1.1185 - BinHeap<Value, EdgeIntMap> *_delta3;
1.1186 + IntEdgeMap *_delta3_index;
1.1187 + BinHeap<Value, IntEdgeMap> *_delta3;
1.1188
1.1189 IntIntMap *_delta4_index;
1.1190 BinHeap<Value, IntIntMap> *_delta4;
1.1191 @@ -775,14 +765,14 @@
1.1192 _node_potential = new NodePotential(_graph);
1.1193 }
1.1194 if (!_blossom_set) {
1.1195 - _blossom_index = new NodeIntMap(_graph);
1.1196 + _blossom_index = new IntNodeMap(_graph);
1.1197 _blossom_set = new BlossomSet(*_blossom_index);
1.1198 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
1.1199 }
1.1200
1.1201 if (!_node_index) {
1.1202 - _node_index = new NodeIntMap(_graph);
1.1203 - _node_heap_index = new ArcIntMap(_graph);
1.1204 + _node_index = new IntNodeMap(_graph);
1.1205 + _node_heap_index = new IntArcMap(_graph);
1.1206 _node_data = new RangeMap<NodeData>(_node_num,
1.1207 NodeData(*_node_heap_index));
1.1208 }
1.1209 @@ -792,16 +782,16 @@
1.1210 _tree_set = new TreeSet(*_tree_set_index);
1.1211 }
1.1212 if (!_delta1) {
1.1213 - _delta1_index = new NodeIntMap(_graph);
1.1214 - _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);
1.1215 + _delta1_index = new IntNodeMap(_graph);
1.1216 + _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
1.1217 }
1.1218 if (!_delta2) {
1.1219 _delta2_index = new IntIntMap(_blossom_num);
1.1220 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
1.1221 }
1.1222 if (!_delta3) {
1.1223 - _delta3_index = new EdgeIntMap(_graph);
1.1224 - _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
1.1225 + _delta3_index = new IntEdgeMap(_graph);
1.1226 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1.1227 }
1.1228 if (!_delta4) {
1.1229 _delta4_index = new IntIntMap(_blossom_num);
1.1230 @@ -1266,10 +1256,10 @@
1.1231 }
1.1232
1.1233
1.1234 - void augmentOnArc(const Edge& arc) {
1.1235 -
1.1236 - int left = _blossom_set->find(_graph.u(arc));
1.1237 - int right = _blossom_set->find(_graph.v(arc));
1.1238 + void augmentOnEdge(const Edge& edge) {
1.1239 +
1.1240 + int left = _blossom_set->find(_graph.u(edge));
1.1241 + int right = _blossom_set->find(_graph.v(edge));
1.1242
1.1243 if ((*_blossom_data)[left].status == EVEN) {
1.1244 int left_tree = _tree_set->find(left);
1.1245 @@ -1289,8 +1279,8 @@
1.1246 unmatchedToMatched(right);
1.1247 }
1.1248
1.1249 - (*_blossom_data)[left].next = _graph.direct(arc, true);
1.1250 - (*_blossom_data)[right].next = _graph.direct(arc, false);
1.1251 + (*_blossom_data)[left].next = _graph.direct(edge, true);
1.1252 + (*_blossom_data)[right].next = _graph.direct(edge, false);
1.1253 }
1.1254
1.1255 void extendOnArc(const Arc& arc) {
1.1256 @@ -1310,7 +1300,7 @@
1.1257 matchedToEven(even, tree);
1.1258 }
1.1259
1.1260 - void shrinkOnArc(const Edge& edge, int tree) {
1.1261 + void shrinkOnEdge(const Edge& edge, int tree) {
1.1262 int nca = -1;
1.1263 std::vector<int> left_path, right_path;
1.1264
1.1265 @@ -1652,7 +1642,7 @@
1.1266 createStructures();
1.1267
1.1268 for (ArcIt e(_graph); e != INVALID; ++e) {
1.1269 - _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
1.1270 + _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
1.1271 }
1.1272 for (NodeIt n(_graph); n != INVALID; ++n) {
1.1273 _delta1_index->set(n, _delta1->PRE_HEAP);
1.1274 @@ -1769,9 +1759,9 @@
1.1275 }
1.1276
1.1277 if (left_tree == right_tree) {
1.1278 - shrinkOnArc(e, left_tree);
1.1279 + shrinkOnEdge(e, left_tree);
1.1280 } else {
1.1281 - augmentOnArc(e);
1.1282 + augmentOnEdge(e);
1.1283 unmatched -= 2;
1.1284 }
1.1285 }
1.1286 @@ -1818,11 +1808,24 @@
1.1287 return sum /= 2;
1.1288 }
1.1289
1.1290 - /// \brief Returns true when the arc is in the matching.
1.1291 + /// \brief Returns the cardinality of the matching.
1.1292 ///
1.1293 - /// Returns true when the arc is in the matching.
1.1294 - bool matching(const Edge& arc) const {
1.1295 - return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
1.1296 + /// Returns the cardinality of the matching.
1.1297 + int matchingSize() const {
1.1298 + int num = 0;
1.1299 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1300 + if ((*_matching)[n] != INVALID) {
1.1301 + ++num;
1.1302 + }
1.1303 + }
1.1304 + return num /= 2;
1.1305 + }
1.1306 +
1.1307 + /// \brief Returns true when the edge is in the matching.
1.1308 + ///
1.1309 + /// Returns true when the edge is in the matching.
1.1310 + bool matching(const Edge& edge) const {
1.1311 + return edge == (*_matching)[_graph.u(edge)];
1.1312 }
1.1313
1.1314 /// \brief Returns the incident matching arc.
1.1315 @@ -1913,16 +1916,11 @@
1.1316 _last = _algorithm->_blossom_potential[variable].end;
1.1317 }
1.1318
1.1319 - /// \brief Invalid constructor.
1.1320 - ///
1.1321 - /// Invalid constructor.
1.1322 - BlossomIt(Invalid) : _index(-1) {}
1.1323 -
1.1324 /// \brief Conversion to node.
1.1325 ///
1.1326 /// Conversion to node.
1.1327 operator Node() const {
1.1328 - return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
1.1329 + return _algorithm->_blossom_node_list[_index];
1.1330 }
1.1331
1.1332 /// \brief Increment operator.
1.1333 @@ -1930,18 +1928,18 @@
1.1334 /// Increment operator.
1.1335 BlossomIt& operator++() {
1.1336 ++_index;
1.1337 - if (_index == _last) {
1.1338 - _index = -1;
1.1339 - }
1.1340 return *this;
1.1341 }
1.1342
1.1343 - bool operator==(const BlossomIt& it) const {
1.1344 - return _index == it._index;
1.1345 - }
1.1346 - bool operator!=(const BlossomIt& it) const {
1.1347 - return _index != it._index;
1.1348 - }
1.1349 + /// \brief Validity checking
1.1350 + ///
1.1351 + /// Checks whether the iterator is invalid.
1.1352 + bool operator==(Invalid) const { return _index == _last; }
1.1353 +
1.1354 + /// \brief Validity checking
1.1355 + ///
1.1356 + /// Checks whether the iterator is valid.
1.1357 + bool operator!=(Invalid) const { return _index != _last; }
1.1358
1.1359 private:
1.1360 const MaxWeightedMatching* _algorithm;
1.1361 @@ -1958,29 +1956,32 @@
1.1362 /// \brief Weighted perfect matching in general graphs
1.1363 ///
1.1364 /// This class provides an efficient implementation of Edmond's
1.1365 - /// maximum weighted perfecr matching algorithm. The implementation
1.1366 + /// maximum weighted perfect matching algorithm. The implementation
1.1367 /// is based on extensive use of priority queues and provides
1.1368 /// \f$O(nm\log(n))\f$ time complexity.
1.1369 ///
1.1370 /// The maximum weighted matching problem is to find undirected
1.1371 - /// arcs in the digraph with maximum overall weight and no two of
1.1372 - /// them shares their endpoints and covers all nodes. The problem
1.1373 - /// can be formulated with the next linear program:
1.1374 + /// edges in the graph with maximum overall weight and no two of
1.1375 + /// them shares their ends and covers all nodes. The problem can be
1.1376 + /// formulated with the following linear program.
1.1377 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1.1378 - ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
1.1379 + /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
1.1380 + \quad \forall B\in\mathcal{O}\f] */
1.1381 /// \f[x_e \ge 0\quad \forall e\in E\f]
1.1382 /// \f[\max \sum_{e\in E}x_ew_e\f]
1.1383 - /// where \f$\delta(X)\f$ is the set of arcs incident to a node in
1.1384 - /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in
1.1385 - /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
1.1386 - /// the nodes.
1.1387 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1.1388 + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
1.1389 + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
1.1390 + /// subsets of the nodes.
1.1391 ///
1.1392 /// The algorithm calculates an optimal matching and a proof of the
1.1393 /// optimality. The solution of the dual problem can be used to check
1.1394 - /// the result of the algorithm. The dual linear problem is the next:
1.1395 - /// \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]
1.1396 + /// the result of the algorithm. The dual linear problem is the
1.1397 + /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
1.1398 + w_{uv} \quad \forall uv\in E\f] */
1.1399 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1.1400 - /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
1.1401 + /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1.1402 + \frac{\vert B \vert - 1}{2}z_B\f] */
1.1403 ///
1.1404 /// The algorithm can be executed with \c run() or the \c init() and
1.1405 /// then the \c start() member functions. After it the matching can
1.1406 @@ -2040,16 +2041,13 @@
1.1407 int _node_num;
1.1408 int _blossom_num;
1.1409
1.1410 - typedef typename Graph::template NodeMap<int> NodeIntMap;
1.1411 - typedef typename Graph::template ArcMap<int> ArcIntMap;
1.1412 - typedef typename Graph::template EdgeMap<int> EdgeIntMap;
1.1413 typedef RangeMap<int> IntIntMap;
1.1414
1.1415 enum Status {
1.1416 EVEN = -1, MATCHED = 0, ODD = 1
1.1417 };
1.1418
1.1419 - typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
1.1420 + typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
1.1421 struct BlossomData {
1.1422 int tree;
1.1423 Status status;
1.1424 @@ -2057,21 +2055,21 @@
1.1425 Value pot, offset;
1.1426 };
1.1427
1.1428 - NodeIntMap *_blossom_index;
1.1429 + IntNodeMap *_blossom_index;
1.1430 BlossomSet *_blossom_set;
1.1431 RangeMap<BlossomData>* _blossom_data;
1.1432
1.1433 - NodeIntMap *_node_index;
1.1434 - ArcIntMap *_node_heap_index;
1.1435 + IntNodeMap *_node_index;
1.1436 + IntArcMap *_node_heap_index;
1.1437
1.1438 struct NodeData {
1.1439
1.1440 - NodeData(ArcIntMap& node_heap_index)
1.1441 + NodeData(IntArcMap& node_heap_index)
1.1442 : heap(node_heap_index) {}
1.1443
1.1444 int blossom;
1.1445 Value pot;
1.1446 - BinHeap<Value, ArcIntMap> heap;
1.1447 + BinHeap<Value, IntArcMap> heap;
1.1448 std::map<int, Arc> heap_index;
1.1449
1.1450 int tree;
1.1451 @@ -2087,8 +2085,8 @@
1.1452 IntIntMap *_delta2_index;
1.1453 BinHeap<Value, IntIntMap> *_delta2;
1.1454
1.1455 - EdgeIntMap *_delta3_index;
1.1456 - BinHeap<Value, EdgeIntMap> *_delta3;
1.1457 + IntEdgeMap *_delta3_index;
1.1458 + BinHeap<Value, IntEdgeMap> *_delta3;
1.1459
1.1460 IntIntMap *_delta4_index;
1.1461 BinHeap<Value, IntIntMap> *_delta4;
1.1462 @@ -2106,16 +2104,16 @@
1.1463 _node_potential = new NodePotential(_graph);
1.1464 }
1.1465 if (!_blossom_set) {
1.1466 - _blossom_index = new NodeIntMap(_graph);
1.1467 + _blossom_index = new IntNodeMap(_graph);
1.1468 _blossom_set = new BlossomSet(*_blossom_index);
1.1469 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
1.1470 }
1.1471
1.1472 if (!_node_index) {
1.1473 - _node_index = new NodeIntMap(_graph);
1.1474 - _node_heap_index = new ArcIntMap(_graph);
1.1475 + _node_index = new IntNodeMap(_graph);
1.1476 + _node_heap_index = new IntArcMap(_graph);
1.1477 _node_data = new RangeMap<NodeData>(_node_num,
1.1478 - NodeData(*_node_heap_index));
1.1479 + NodeData(*_node_heap_index));
1.1480 }
1.1481
1.1482 if (!_tree_set) {
1.1483 @@ -2127,8 +2125,8 @@
1.1484 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
1.1485 }
1.1486 if (!_delta3) {
1.1487 - _delta3_index = new EdgeIntMap(_graph);
1.1488 - _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
1.1489 + _delta3_index = new IntEdgeMap(_graph);
1.1490 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1.1491 }
1.1492 if (!_delta4) {
1.1493 _delta4_index = new IntIntMap(_blossom_num);
1.1494 @@ -2461,10 +2459,10 @@
1.1495 _tree_set->eraseClass(tree);
1.1496 }
1.1497
1.1498 - void augmentOnArc(const Edge& arc) {
1.1499 -
1.1500 - int left = _blossom_set->find(_graph.u(arc));
1.1501 - int right = _blossom_set->find(_graph.v(arc));
1.1502 + void augmentOnEdge(const Edge& edge) {
1.1503 +
1.1504 + int left = _blossom_set->find(_graph.u(edge));
1.1505 + int right = _blossom_set->find(_graph.v(edge));
1.1506
1.1507 int left_tree = _tree_set->find(left);
1.1508 alternatePath(left, left_tree);
1.1509 @@ -2474,8 +2472,8 @@
1.1510 alternatePath(right, right_tree);
1.1511 destroyTree(right_tree);
1.1512
1.1513 - (*_blossom_data)[left].next = _graph.direct(arc, true);
1.1514 - (*_blossom_data)[right].next = _graph.direct(arc, false);
1.1515 + (*_blossom_data)[left].next = _graph.direct(edge, true);
1.1516 + (*_blossom_data)[right].next = _graph.direct(edge, false);
1.1517 }
1.1518
1.1519 void extendOnArc(const Arc& arc) {
1.1520 @@ -2495,7 +2493,7 @@
1.1521 matchedToEven(even, tree);
1.1522 }
1.1523
1.1524 - void shrinkOnArc(const Edge& edge, int tree) {
1.1525 + void shrinkOnEdge(const Edge& edge, int tree) {
1.1526 int nca = -1;
1.1527 std::vector<int> left_path, right_path;
1.1528
1.1529 @@ -2831,7 +2829,7 @@
1.1530 createStructures();
1.1531
1.1532 for (ArcIt e(_graph); e != INVALID; ++e) {
1.1533 - _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
1.1534 + _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
1.1535 }
1.1536 for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1537 _delta3_index->set(e, _delta3->PRE_HEAP);
1.1538 @@ -2924,9 +2922,9 @@
1.1539 int right_tree = _tree_set->find(right_blossom);
1.1540
1.1541 if (left_tree == right_tree) {
1.1542 - shrinkOnArc(e, left_tree);
1.1543 + shrinkOnEdge(e, left_tree);
1.1544 } else {
1.1545 - augmentOnArc(e);
1.1546 + augmentOnEdge(e);
1.1547 unmatched -= 2;
1.1548 }
1.1549 }
1.1550 @@ -2974,16 +2972,16 @@
1.1551 return sum /= 2;
1.1552 }
1.1553
1.1554 - /// \brief Returns true when the arc is in the matching.
1.1555 + /// \brief Returns true when the edge is in the matching.
1.1556 ///
1.1557 - /// Returns true when the arc is in the matching.
1.1558 - bool matching(const Edge& arc) const {
1.1559 - return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
1.1560 + /// Returns true when the edge is in the matching.
1.1561 + bool matching(const Edge& edge) const {
1.1562 + return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
1.1563 }
1.1564
1.1565 - /// \brief Returns the incident matching arc.
1.1566 + /// \brief Returns the incident matching edge.
1.1567 ///
1.1568 - /// Returns the incident matching arc from given node.
1.1569 + /// Returns the incident matching arc from given edge.
1.1570 Arc matching(const Node& node) const {
1.1571 return (*_matching)[node];
1.1572 }
1.1573 @@ -3066,16 +3064,11 @@
1.1574 _last = _algorithm->_blossom_potential[variable].end;
1.1575 }
1.1576
1.1577 - /// \brief Invalid constructor.
1.1578 - ///
1.1579 - /// Invalid constructor.
1.1580 - BlossomIt(Invalid) : _index(-1) {}
1.1581 -
1.1582 /// \brief Conversion to node.
1.1583 ///
1.1584 /// Conversion to node.
1.1585 operator Node() const {
1.1586 - return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
1.1587 + return _algorithm->_blossom_node_list[_index];
1.1588 }
1.1589
1.1590 /// \brief Increment operator.
1.1591 @@ -3083,18 +3076,18 @@
1.1592 /// Increment operator.
1.1593 BlossomIt& operator++() {
1.1594 ++_index;
1.1595 - if (_index == _last) {
1.1596 - _index = -1;
1.1597 - }
1.1598 return *this;
1.1599 }
1.1600
1.1601 - bool operator==(const BlossomIt& it) const {
1.1602 - return _index == it._index;
1.1603 - }
1.1604 - bool operator!=(const BlossomIt& it) const {
1.1605 - return _index != it._index;
1.1606 - }
1.1607 + /// \brief Validity checking
1.1608 + ///
1.1609 + /// Checks whether the iterator is invalid.
1.1610 + bool operator==(Invalid) const { return _index == _last; }
1.1611 +
1.1612 + /// \brief Validity checking
1.1613 + ///
1.1614 + /// Checks whether the iterator is valid.
1.1615 + bool operator!=(Invalid) const { return _index != _last; }
1.1616
1.1617 private:
1.1618 const MaxWeightedPerfectMatching* _algorithm;