COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/matching.h @ 642:111698359429

Last change on this file since 642:111698359429 was 594:d657c71db7db, checked in by Peter Kovacs <kpeter@…>, 15 years ago

Rename max_matching.h to matching.h (#265)

File size: 102.2 KB
Line 
1/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library.
4 *
5 * Copyright (C) 2003-2009
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_MAX_MATCHING_H
20#define LEMON_MAX_MATCHING_H
21
22#include <vector>
23#include <queue>
24#include <set>
25#include <limits>
26
27#include <lemon/core.h>
28#include <lemon/unionfind.h>
29#include <lemon/bin_heap.h>
30#include <lemon/maps.h>
31
32///\ingroup matching
33///\file
34///\brief Maximum matching algorithms in general graphs.
35
36namespace lemon {
37
38  /// \ingroup matching
39  ///
40  /// \brief Maximum cardinality matching in general graphs
41  ///
42  /// This class implements Edmonds' alternating forest matching algorithm
43  /// for finding a maximum cardinality matching in a general undirected graph.
44  /// It can be started from an arbitrary initial matching
45  /// (the default is the empty one).
46  ///
47  /// The dual solution of the problem is a map of the nodes to
48  /// \ref MaxMatching::Status "Status", having values \c EVEN (or \c D),
49  /// \c ODD (or \c A) and \c MATCHED (or \c C) defining the Gallai-Edmonds
50  /// decomposition of the graph. The nodes in \c EVEN/D induce a subgraph
51  /// with factor-critical components, the nodes in \c ODD/A form the
52  /// canonical barrier, and the nodes in \c MATCHED/C induce a graph having
53  /// a perfect matching. The number of the factor-critical components
54  /// minus the number of barrier nodes is a lower bound on the
55  /// unmatched nodes, and the matching is optimal if and only if this bound is
56  /// tight. This decomposition can be obtained using \ref status() or
57  /// \ref statusMap() after running the algorithm.
58  ///
59  /// \tparam GR The undirected graph type the algorithm runs on.
60  template <typename GR>
61  class MaxMatching {
62  public:
63
64    /// The graph type of the algorithm
65    typedef GR Graph;
66    /// The type of the matching map
67    typedef typename Graph::template NodeMap<typename Graph::Arc>
68    MatchingMap;
69
70    ///\brief Status constants for Gallai-Edmonds decomposition.
71    ///
72    ///These constants are used for indicating the Gallai-Edmonds
73    ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
74    ///induce a subgraph with factor-critical components, the nodes with
75    ///status \c ODD (or \c A) form the canonical barrier, and the nodes
76    ///with status \c MATCHED (or \c C) induce a subgraph having a
77    ///perfect matching.
78    enum Status {
79      EVEN = 1,       ///< = 1. (\c D is an alias for \c EVEN.)
80      D = 1,
81      MATCHED = 0,    ///< = 0. (\c C is an alias for \c MATCHED.)
82      C = 0,
83      ODD = -1,       ///< = -1. (\c A is an alias for \c ODD.)
84      A = -1,
85      UNMATCHED = -2  ///< = -2.
86    };
87
88    /// The type of the status map
89    typedef typename Graph::template NodeMap<Status> StatusMap;
90
91  private:
92
93    TEMPLATE_GRAPH_TYPEDEFS(Graph);
94
95    typedef UnionFindEnum<IntNodeMap> BlossomSet;
96    typedef ExtendFindEnum<IntNodeMap> TreeSet;
97    typedef RangeMap<Node> NodeIntMap;
98    typedef MatchingMap EarMap;
99    typedef std::vector<Node> NodeQueue;
100
101    const Graph& _graph;
102    MatchingMap* _matching;
103    StatusMap* _status;
104
105    EarMap* _ear;
106
107    IntNodeMap* _blossom_set_index;
108    BlossomSet* _blossom_set;
109    NodeIntMap* _blossom_rep;
110
111    IntNodeMap* _tree_set_index;
112    TreeSet* _tree_set;
113
114    NodeQueue _node_queue;
115    int _process, _postpone, _last;
116
117    int _node_num;
118
119  private:
120
121    void createStructures() {
122      _node_num = countNodes(_graph);
123      if (!_matching) {
124        _matching = new MatchingMap(_graph);
125      }
126      if (!_status) {
127        _status = new StatusMap(_graph);
128      }
129      if (!_ear) {
130        _ear = new EarMap(_graph);
131      }
132      if (!_blossom_set) {
133        _blossom_set_index = new IntNodeMap(_graph);
134        _blossom_set = new BlossomSet(*_blossom_set_index);
135      }
136      if (!_blossom_rep) {
137        _blossom_rep = new NodeIntMap(_node_num);
138      }
139      if (!_tree_set) {
140        _tree_set_index = new IntNodeMap(_graph);
141        _tree_set = new TreeSet(*_tree_set_index);
142      }
143      _node_queue.resize(_node_num);
144    }
145
146    void destroyStructures() {
147      if (_matching) {
148        delete _matching;
149      }
150      if (_status) {
151        delete _status;
152      }
153      if (_ear) {
154        delete _ear;
155      }
156      if (_blossom_set) {
157        delete _blossom_set;
158        delete _blossom_set_index;
159      }
160      if (_blossom_rep) {
161        delete _blossom_rep;
162      }
163      if (_tree_set) {
164        delete _tree_set_index;
165        delete _tree_set;
166      }
167    }
168
169    void processDense(const Node& n) {
170      _process = _postpone = _last = 0;
171      _node_queue[_last++] = n;
172
173      while (_process != _last) {
174        Node u = _node_queue[_process++];
175        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
176          Node v = _graph.target(a);
177          if ((*_status)[v] == MATCHED) {
178            extendOnArc(a);
179          } else if ((*_status)[v] == UNMATCHED) {
180            augmentOnArc(a);
181            return;
182          }
183        }
184      }
185
186      while (_postpone != _last) {
187        Node u = _node_queue[_postpone++];
188
189        for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
190          Node v = _graph.target(a);
191
192          if ((*_status)[v] == EVEN) {
193            if (_blossom_set->find(u) != _blossom_set->find(v)) {
194              shrinkOnEdge(a);
195            }
196          }
197
198          while (_process != _last) {
199            Node w = _node_queue[_process++];
200            for (OutArcIt b(_graph, w); b != INVALID; ++b) {
201              Node x = _graph.target(b);
202              if ((*_status)[x] == MATCHED) {
203                extendOnArc(b);
204              } else if ((*_status)[x] == UNMATCHED) {
205                augmentOnArc(b);
206                return;
207              }
208            }
209          }
210        }
211      }
212    }
213
214    void processSparse(const Node& n) {
215      _process = _last = 0;
216      _node_queue[_last++] = n;
217      while (_process != _last) {
218        Node u = _node_queue[_process++];
219        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
220          Node v = _graph.target(a);
221
222          if ((*_status)[v] == EVEN) {
223            if (_blossom_set->find(u) != _blossom_set->find(v)) {
224              shrinkOnEdge(a);
225            }
226          } else if ((*_status)[v] == MATCHED) {
227            extendOnArc(a);
228          } else if ((*_status)[v] == UNMATCHED) {
229            augmentOnArc(a);
230            return;
231          }
232        }
233      }
234    }
235
236    void shrinkOnEdge(const Edge& e) {
237      Node nca = INVALID;
238
239      {
240        std::set<Node> left_set, right_set;
241
242        Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
243        left_set.insert(left);
244
245        Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
246        right_set.insert(right);
247
248        while (true) {
249          if ((*_matching)[left] == INVALID) break;
250          left = _graph.target((*_matching)[left]);
251          left = (*_blossom_rep)[_blossom_set->
252                                 find(_graph.target((*_ear)[left]))];
253          if (right_set.find(left) != right_set.end()) {
254            nca = left;
255            break;
256          }
257          left_set.insert(left);
258
259          if ((*_matching)[right] == INVALID) break;
260          right = _graph.target((*_matching)[right]);
261          right = (*_blossom_rep)[_blossom_set->
262                                  find(_graph.target((*_ear)[right]))];
263          if (left_set.find(right) != left_set.end()) {
264            nca = right;
265            break;
266          }
267          right_set.insert(right);
268        }
269
270        if (nca == INVALID) {
271          if ((*_matching)[left] == INVALID) {
272            nca = right;
273            while (left_set.find(nca) == left_set.end()) {
274              nca = _graph.target((*_matching)[nca]);
275              nca =(*_blossom_rep)[_blossom_set->
276                                   find(_graph.target((*_ear)[nca]))];
277            }
278          } else {
279            nca = left;
280            while (right_set.find(nca) == right_set.end()) {
281              nca = _graph.target((*_matching)[nca]);
282              nca = (*_blossom_rep)[_blossom_set->
283                                   find(_graph.target((*_ear)[nca]))];
284            }
285          }
286        }
287      }
288
289      {
290
291        Node node = _graph.u(e);
292        Arc arc = _graph.direct(e, true);
293        Node base = (*_blossom_rep)[_blossom_set->find(node)];
294
295        while (base != nca) {
296          (*_ear)[node] = arc;
297
298          Node n = node;
299          while (n != base) {
300            n = _graph.target((*_matching)[n]);
301            Arc a = (*_ear)[n];
302            n = _graph.target(a);
303            (*_ear)[n] = _graph.oppositeArc(a);
304          }
305          node = _graph.target((*_matching)[base]);
306          _tree_set->erase(base);
307          _tree_set->erase(node);
308          _blossom_set->insert(node, _blossom_set->find(base));
309          (*_status)[node] = EVEN;
310          _node_queue[_last++] = node;
311          arc = _graph.oppositeArc((*_ear)[node]);
312          node = _graph.target((*_ear)[node]);
313          base = (*_blossom_rep)[_blossom_set->find(node)];
314          _blossom_set->join(_graph.target(arc), base);
315        }
316      }
317
318      (*_blossom_rep)[_blossom_set->find(nca)] = nca;
319
320      {
321
322        Node node = _graph.v(e);
323        Arc arc = _graph.direct(e, false);
324        Node base = (*_blossom_rep)[_blossom_set->find(node)];
325
326        while (base != nca) {
327          (*_ear)[node] = arc;
328
329          Node n = node;
330          while (n != base) {
331            n = _graph.target((*_matching)[n]);
332            Arc a = (*_ear)[n];
333            n = _graph.target(a);
334            (*_ear)[n] = _graph.oppositeArc(a);
335          }
336          node = _graph.target((*_matching)[base]);
337          _tree_set->erase(base);
338          _tree_set->erase(node);
339          _blossom_set->insert(node, _blossom_set->find(base));
340          (*_status)[node] = EVEN;
341          _node_queue[_last++] = node;
342          arc = _graph.oppositeArc((*_ear)[node]);
343          node = _graph.target((*_ear)[node]);
344          base = (*_blossom_rep)[_blossom_set->find(node)];
345          _blossom_set->join(_graph.target(arc), base);
346        }
347      }
348
349      (*_blossom_rep)[_blossom_set->find(nca)] = nca;
350    }
351
352    void extendOnArc(const Arc& a) {
353      Node base = _graph.source(a);
354      Node odd = _graph.target(a);
355
356      (*_ear)[odd] = _graph.oppositeArc(a);
357      Node even = _graph.target((*_matching)[odd]);
358      (*_blossom_rep)[_blossom_set->insert(even)] = even;
359      (*_status)[odd] = ODD;
360      (*_status)[even] = EVEN;
361      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
362      _tree_set->insert(odd, tree);
363      _tree_set->insert(even, tree);
364      _node_queue[_last++] = even;
365
366    }
367
368    void augmentOnArc(const Arc& a) {
369      Node even = _graph.source(a);
370      Node odd = _graph.target(a);
371
372      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
373
374      (*_matching)[odd] = _graph.oppositeArc(a);
375      (*_status)[odd] = MATCHED;
376
377      Arc arc = (*_matching)[even];
378      (*_matching)[even] = a;
379
380      while (arc != INVALID) {
381        odd = _graph.target(arc);
382        arc = (*_ear)[odd];
383        even = _graph.target(arc);
384        (*_matching)[odd] = arc;
385        arc = (*_matching)[even];
386        (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]);
387      }
388
389      for (typename TreeSet::ItemIt it(*_tree_set, tree);
390           it != INVALID; ++it) {
391        if ((*_status)[it] == ODD) {
392          (*_status)[it] = MATCHED;
393        } else {
394          int blossom = _blossom_set->find(it);
395          for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
396               jt != INVALID; ++jt) {
397            (*_status)[jt] = MATCHED;
398          }
399          _blossom_set->eraseClass(blossom);
400        }
401      }
402      _tree_set->eraseClass(tree);
403
404    }
405
406  public:
407
408    /// \brief Constructor
409    ///
410    /// Constructor.
411    MaxMatching(const Graph& graph)
412      : _graph(graph), _matching(0), _status(0), _ear(0),
413        _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
414        _tree_set_index(0), _tree_set(0) {}
415
416    ~MaxMatching() {
417      destroyStructures();
418    }
419
420    /// \name Execution Control
421    /// The simplest way to execute the algorithm is to use the
422    /// \c run() member function.\n
423    /// If you need better control on the execution, you have to call
424    /// one of the functions \ref init(), \ref greedyInit() or
425    /// \ref matchingInit() first, then you can start the algorithm with
426    /// \ref startSparse() or \ref startDense().
427
428    ///@{
429
430    /// \brief Set the initial matching to the empty matching.
431    ///
432    /// This function sets the initial matching to the empty matching.
433    void init() {
434      createStructures();
435      for(NodeIt n(_graph); n != INVALID; ++n) {
436        (*_matching)[n] = INVALID;
437        (*_status)[n] = UNMATCHED;
438      }
439    }
440
441    /// \brief Find an initial matching in a greedy way.
442    ///
443    /// This function finds an initial matching in a greedy way.
444    void greedyInit() {
445      createStructures();
446      for (NodeIt n(_graph); n != INVALID; ++n) {
447        (*_matching)[n] = INVALID;
448        (*_status)[n] = UNMATCHED;
449      }
450      for (NodeIt n(_graph); n != INVALID; ++n) {
451        if ((*_matching)[n] == INVALID) {
452          for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
453            Node v = _graph.target(a);
454            if ((*_matching)[v] == INVALID && v != n) {
455              (*_matching)[n] = a;
456              (*_status)[n] = MATCHED;
457              (*_matching)[v] = _graph.oppositeArc(a);
458              (*_status)[v] = MATCHED;
459              break;
460            }
461          }
462        }
463      }
464    }
465
466
467    /// \brief Initialize the matching from a map.
468    ///
469    /// This function initializes the matching from a \c bool valued edge
470    /// map. This map should have the property that there are no two incident
471    /// edges with \c true value, i.e. it really contains a matching.
472    /// \return \c true if the map contains a matching.
473    template <typename MatchingMap>
474    bool matchingInit(const MatchingMap& matching) {
475      createStructures();
476
477      for (NodeIt n(_graph); n != INVALID; ++n) {
478        (*_matching)[n] = INVALID;
479        (*_status)[n] = UNMATCHED;
480      }
481      for(EdgeIt e(_graph); e!=INVALID; ++e) {
482        if (matching[e]) {
483
484          Node u = _graph.u(e);
485          if ((*_matching)[u] != INVALID) return false;
486          (*_matching)[u] = _graph.direct(e, true);
487          (*_status)[u] = MATCHED;
488
489          Node v = _graph.v(e);
490          if ((*_matching)[v] != INVALID) return false;
491          (*_matching)[v] = _graph.direct(e, false);
492          (*_status)[v] = MATCHED;
493        }
494      }
495      return true;
496    }
497
498    /// \brief Start Edmonds' algorithm
499    ///
500    /// This function runs the original Edmonds' algorithm.
501    ///
502    /// \pre \ref Init(), \ref greedyInit() or \ref matchingInit() must be
503    /// called before using this function.
504    void startSparse() {
505      for(NodeIt n(_graph); n != INVALID; ++n) {
506        if ((*_status)[n] == UNMATCHED) {
507          (*_blossom_rep)[_blossom_set->insert(n)] = n;
508          _tree_set->insert(n);
509          (*_status)[n] = EVEN;
510          processSparse(n);
511        }
512      }
513    }
514
515    /// \brief Start Edmonds' algorithm with a heuristic improvement
516    /// for dense graphs
517    ///
518    /// This function runs Edmonds' algorithm with a heuristic of postponing
519    /// shrinks, therefore resulting in a faster algorithm for dense graphs.
520    ///
521    /// \pre \ref Init(), \ref greedyInit() or \ref matchingInit() must be
522    /// called before using this function.
523    void startDense() {
524      for(NodeIt n(_graph); n != INVALID; ++n) {
525        if ((*_status)[n] == UNMATCHED) {
526          (*_blossom_rep)[_blossom_set->insert(n)] = n;
527          _tree_set->insert(n);
528          (*_status)[n] = EVEN;
529          processDense(n);
530        }
531      }
532    }
533
534
535    /// \brief Run Edmonds' algorithm
536    ///
537    /// This function runs Edmonds' algorithm. An additional heuristic of
538    /// postponing shrinks is used for relatively dense graphs
539    /// (for which <tt>m>=2*n</tt> holds).
540    void run() {
541      if (countEdges(_graph) < 2 * countNodes(_graph)) {
542        greedyInit();
543        startSparse();
544      } else {
545        init();
546        startDense();
547      }
548    }
549
550    /// @}
551
552    /// \name Primal Solution
553    /// Functions to get the primal solution, i.e. the maximum matching.
554
555    /// @{
556
557    /// \brief Return the size (cardinality) of the matching.
558    ///
559    /// This function returns the size (cardinality) of the current matching.
560    /// After run() it returns the size of the maximum matching in the graph.
561    int matchingSize() const {
562      int size = 0;
563      for (NodeIt n(_graph); n != INVALID; ++n) {
564        if ((*_matching)[n] != INVALID) {
565          ++size;
566        }
567      }
568      return size / 2;
569    }
570
571    /// \brief Return \c true if the given edge is in the matching.
572    ///
573    /// This function returns \c true if the given edge is in the current
574    /// matching.
575    bool matching(const Edge& edge) const {
576      return edge == (*_matching)[_graph.u(edge)];
577    }
578
579    /// \brief Return the matching arc (or edge) incident to the given node.
580    ///
581    /// This function returns the matching arc (or edge) incident to the
582    /// given node in the current matching or \c INVALID if the node is
583    /// not covered by the matching.
584    Arc matching(const Node& n) const {
585      return (*_matching)[n];
586    }
587
588    /// \brief Return a const reference to the matching map.
589    ///
590    /// This function returns a const reference to a node map that stores
591    /// the matching arc (or edge) incident to each node.
592    const MatchingMap& matchingMap() const {
593      return *_matching;
594    }
595
596    /// \brief Return the mate of the given node.
597    ///
598    /// This function returns the mate of the given node in the current
599    /// matching or \c INVALID if the node is not covered by the matching.
600    Node mate(const Node& n) const {
601      return (*_matching)[n] != INVALID ?
602        _graph.target((*_matching)[n]) : INVALID;
603    }
604
605    /// @}
606
607    /// \name Dual Solution
608    /// Functions to get the dual solution, i.e. the Gallai-Edmonds
609    /// decomposition.
610
611    /// @{
612
613    /// \brief Return the status of the given node in the Edmonds-Gallai
614    /// decomposition.
615    ///
616    /// This function returns the \ref Status "status" of the given node
617    /// in the Edmonds-Gallai decomposition.
618    Status status(const Node& n) const {
619      return (*_status)[n];
620    }
621
622    /// \brief Return a const reference to the status map, which stores
623    /// the Edmonds-Gallai decomposition.
624    ///
625    /// This function returns a const reference to a node map that stores the
626    /// \ref Status "status" of each node in the Edmonds-Gallai decomposition.
627    const StatusMap& statusMap() const {
628      return *_status;
629    }
630
631    /// \brief Return \c true if the given node is in the barrier.
632    ///
633    /// This function returns \c true if the given node is in the barrier.
634    bool barrier(const Node& n) const {
635      return (*_status)[n] == ODD;
636    }
637
638    /// @}
639
640  };
641
642  /// \ingroup matching
643  ///
644  /// \brief Weighted matching in general graphs
645  ///
646  /// This class provides an efficient implementation of Edmond's
647  /// maximum weighted matching algorithm. The implementation is based
648  /// on extensive use of priority queues and provides
649  /// \f$O(nm\log n)\f$ time complexity.
650  ///
651  /// The maximum weighted matching problem is to find a subset of the
652  /// edges in an undirected graph with maximum overall weight for which
653  /// each node has at most one incident edge.
654  /// It can be formulated with the following linear program.
655  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
656  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
657      \quad \forall B\in\mathcal{O}\f] */
658  /// \f[x_e \ge 0\quad \forall e\in E\f]
659  /// \f[\max \sum_{e\in E}x_ew_e\f]
660  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
661  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
662  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
663  /// subsets of the nodes.
664  ///
665  /// The algorithm calculates an optimal matching and a proof of the
666  /// optimality. The solution of the dual problem can be used to check
667  /// the result of the algorithm. The dual linear problem is the
668  /// following.
669  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
670      z_B \ge w_{uv} \quad \forall uv\in E\f] */
671  /// \f[y_u \ge 0 \quad \forall u \in V\f]
672  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
673  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
674      \frac{\vert B \vert - 1}{2}z_B\f] */
675  ///
676  /// The algorithm can be executed with the run() function.
677  /// After it the matching (the primal solution) and the dual solution
678  /// can be obtained using the query functions and the
679  /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
680  /// which is able to iterate on the nodes of a blossom.
681  /// If the value type is integer, then the dual solution is multiplied
682  /// by \ref MaxWeightedMatching::dualScale "4".
683  ///
684  /// \tparam GR The undirected graph type the algorithm runs on.
685  /// \tparam WM The type edge weight map. The default type is
686  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
687#ifdef DOXYGEN
688  template <typename GR, typename WM>
689#else
690  template <typename GR,
691            typename WM = typename GR::template EdgeMap<int> >
692#endif
693  class MaxWeightedMatching {
694  public:
695
696    /// The graph type of the algorithm
697    typedef GR Graph;
698    /// The type of the edge weight map
699    typedef WM WeightMap;
700    /// The value type of the edge weights
701    typedef typename WeightMap::Value Value;
702
703    /// The type of the matching map
704    typedef typename Graph::template NodeMap<typename Graph::Arc>
705    MatchingMap;
706
707    /// \brief Scaling factor for dual solution
708    ///
709    /// Scaling factor for dual solution. It is equal to 4 or 1
710    /// according to the value type.
711    static const int dualScale =
712      std::numeric_limits<Value>::is_integer ? 4 : 1;
713
714  private:
715
716    TEMPLATE_GRAPH_TYPEDEFS(Graph);
717
718    typedef typename Graph::template NodeMap<Value> NodePotential;
719    typedef std::vector<Node> BlossomNodeList;
720
721    struct BlossomVariable {
722      int begin, end;
723      Value value;
724
725      BlossomVariable(int _begin, int _end, Value _value)
726        : begin(_begin), end(_end), value(_value) {}
727
728    };
729
730    typedef std::vector<BlossomVariable> BlossomPotential;
731
732    const Graph& _graph;
733    const WeightMap& _weight;
734
735    MatchingMap* _matching;
736
737    NodePotential* _node_potential;
738
739    BlossomPotential _blossom_potential;
740    BlossomNodeList _blossom_node_list;
741
742    int _node_num;
743    int _blossom_num;
744
745    typedef RangeMap<int> IntIntMap;
746
747    enum Status {
748      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
749    };
750
751    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
752    struct BlossomData {
753      int tree;
754      Status status;
755      Arc pred, next;
756      Value pot, offset;
757      Node base;
758    };
759
760    IntNodeMap *_blossom_index;
761    BlossomSet *_blossom_set;
762    RangeMap<BlossomData>* _blossom_data;
763
764    IntNodeMap *_node_index;
765    IntArcMap *_node_heap_index;
766
767    struct NodeData {
768
769      NodeData(IntArcMap& node_heap_index)
770        : heap(node_heap_index) {}
771
772      int blossom;
773      Value pot;
774      BinHeap<Value, IntArcMap> heap;
775      std::map<int, Arc> heap_index;
776
777      int tree;
778    };
779
780    RangeMap<NodeData>* _node_data;
781
782    typedef ExtendFindEnum<IntIntMap> TreeSet;
783
784    IntIntMap *_tree_set_index;
785    TreeSet *_tree_set;
786
787    IntNodeMap *_delta1_index;
788    BinHeap<Value, IntNodeMap> *_delta1;
789
790    IntIntMap *_delta2_index;
791    BinHeap<Value, IntIntMap> *_delta2;
792
793    IntEdgeMap *_delta3_index;
794    BinHeap<Value, IntEdgeMap> *_delta3;
795
796    IntIntMap *_delta4_index;
797    BinHeap<Value, IntIntMap> *_delta4;
798
799    Value _delta_sum;
800
801    void createStructures() {
802      _node_num = countNodes(_graph);
803      _blossom_num = _node_num * 3 / 2;
804
805      if (!_matching) {
806        _matching = new MatchingMap(_graph);
807      }
808      if (!_node_potential) {
809        _node_potential = new NodePotential(_graph);
810      }
811      if (!_blossom_set) {
812        _blossom_index = new IntNodeMap(_graph);
813        _blossom_set = new BlossomSet(*_blossom_index);
814        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
815      }
816
817      if (!_node_index) {
818        _node_index = new IntNodeMap(_graph);
819        _node_heap_index = new IntArcMap(_graph);
820        _node_data = new RangeMap<NodeData>(_node_num,
821                                              NodeData(*_node_heap_index));
822      }
823
824      if (!_tree_set) {
825        _tree_set_index = new IntIntMap(_blossom_num);
826        _tree_set = new TreeSet(*_tree_set_index);
827      }
828      if (!_delta1) {
829        _delta1_index = new IntNodeMap(_graph);
830        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
831      }
832      if (!_delta2) {
833        _delta2_index = new IntIntMap(_blossom_num);
834        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
835      }
836      if (!_delta3) {
837        _delta3_index = new IntEdgeMap(_graph);
838        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
839      }
840      if (!_delta4) {
841        _delta4_index = new IntIntMap(_blossom_num);
842        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
843      }
844    }
845
846    void destroyStructures() {
847      _node_num = countNodes(_graph);
848      _blossom_num = _node_num * 3 / 2;
849
850      if (_matching) {
851        delete _matching;
852      }
853      if (_node_potential) {
854        delete _node_potential;
855      }
856      if (_blossom_set) {
857        delete _blossom_index;
858        delete _blossom_set;
859        delete _blossom_data;
860      }
861
862      if (_node_index) {
863        delete _node_index;
864        delete _node_heap_index;
865        delete _node_data;
866      }
867
868      if (_tree_set) {
869        delete _tree_set_index;
870        delete _tree_set;
871      }
872      if (_delta1) {
873        delete _delta1_index;
874        delete _delta1;
875      }
876      if (_delta2) {
877        delete _delta2_index;
878        delete _delta2;
879      }
880      if (_delta3) {
881        delete _delta3_index;
882        delete _delta3;
883      }
884      if (_delta4) {
885        delete _delta4_index;
886        delete _delta4;
887      }
888    }
889
890    void matchedToEven(int blossom, int tree) {
891      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
892        _delta2->erase(blossom);
893      }
894
895      if (!_blossom_set->trivial(blossom)) {
896        (*_blossom_data)[blossom].pot -=
897          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
898      }
899
900      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
901           n != INVALID; ++n) {
902
903        _blossom_set->increase(n, std::numeric_limits<Value>::max());
904        int ni = (*_node_index)[n];
905
906        (*_node_data)[ni].heap.clear();
907        (*_node_data)[ni].heap_index.clear();
908
909        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
910
911        _delta1->push(n, (*_node_data)[ni].pot);
912
913        for (InArcIt e(_graph, n); e != INVALID; ++e) {
914          Node v = _graph.source(e);
915          int vb = _blossom_set->find(v);
916          int vi = (*_node_index)[v];
917
918          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
919            dualScale * _weight[e];
920
921          if ((*_blossom_data)[vb].status == EVEN) {
922            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
923              _delta3->push(e, rw / 2);
924            }
925          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
926            if (_delta3->state(e) != _delta3->IN_HEAP) {
927              _delta3->push(e, rw);
928            }
929          } else {
930            typename std::map<int, Arc>::iterator it =
931              (*_node_data)[vi].heap_index.find(tree);
932
933            if (it != (*_node_data)[vi].heap_index.end()) {
934              if ((*_node_data)[vi].heap[it->second] > rw) {
935                (*_node_data)[vi].heap.replace(it->second, e);
936                (*_node_data)[vi].heap.decrease(e, rw);
937                it->second = e;
938              }
939            } else {
940              (*_node_data)[vi].heap.push(e, rw);
941              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
942            }
943
944            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
945              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
946
947              if ((*_blossom_data)[vb].status == MATCHED) {
948                if (_delta2->state(vb) != _delta2->IN_HEAP) {
949                  _delta2->push(vb, _blossom_set->classPrio(vb) -
950                               (*_blossom_data)[vb].offset);
951                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
952                           (*_blossom_data)[vb].offset){
953                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
954                                   (*_blossom_data)[vb].offset);
955                }
956              }
957            }
958          }
959        }
960      }
961      (*_blossom_data)[blossom].offset = 0;
962    }
963
964    void matchedToOdd(int blossom) {
965      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
966        _delta2->erase(blossom);
967      }
968      (*_blossom_data)[blossom].offset += _delta_sum;
969      if (!_blossom_set->trivial(blossom)) {
970        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
971                     (*_blossom_data)[blossom].offset);
972      }
973    }
974
975    void evenToMatched(int blossom, int tree) {
976      if (!_blossom_set->trivial(blossom)) {
977        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
978      }
979
980      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
981           n != INVALID; ++n) {
982        int ni = (*_node_index)[n];
983        (*_node_data)[ni].pot -= _delta_sum;
984
985        _delta1->erase(n);
986
987        for (InArcIt e(_graph, n); e != INVALID; ++e) {
988          Node v = _graph.source(e);
989          int vb = _blossom_set->find(v);
990          int vi = (*_node_index)[v];
991
992          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
993            dualScale * _weight[e];
994
995          if (vb == blossom) {
996            if (_delta3->state(e) == _delta3->IN_HEAP) {
997              _delta3->erase(e);
998            }
999          } else if ((*_blossom_data)[vb].status == EVEN) {
1000
1001            if (_delta3->state(e) == _delta3->IN_HEAP) {
1002              _delta3->erase(e);
1003            }
1004
1005            int vt = _tree_set->find(vb);
1006
1007            if (vt != tree) {
1008
1009              Arc r = _graph.oppositeArc(e);
1010
1011              typename std::map<int, Arc>::iterator it =
1012                (*_node_data)[ni].heap_index.find(vt);
1013
1014              if (it != (*_node_data)[ni].heap_index.end()) {
1015                if ((*_node_data)[ni].heap[it->second] > rw) {
1016                  (*_node_data)[ni].heap.replace(it->second, r);
1017                  (*_node_data)[ni].heap.decrease(r, rw);
1018                  it->second = r;
1019                }
1020              } else {
1021                (*_node_data)[ni].heap.push(r, rw);
1022                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1023              }
1024
1025              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1026                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1027
1028                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1029                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1030                               (*_blossom_data)[blossom].offset);
1031                } else if ((*_delta2)[blossom] >
1032                           _blossom_set->classPrio(blossom) -
1033                           (*_blossom_data)[blossom].offset){
1034                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1035                                   (*_blossom_data)[blossom].offset);
1036                }
1037              }
1038            }
1039
1040          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1041            if (_delta3->state(e) == _delta3->IN_HEAP) {
1042              _delta3->erase(e);
1043            }
1044          } else {
1045
1046            typename std::map<int, Arc>::iterator it =
1047              (*_node_data)[vi].heap_index.find(tree);
1048
1049            if (it != (*_node_data)[vi].heap_index.end()) {
1050              (*_node_data)[vi].heap.erase(it->second);
1051              (*_node_data)[vi].heap_index.erase(it);
1052              if ((*_node_data)[vi].heap.empty()) {
1053                _blossom_set->increase(v, std::numeric_limits<Value>::max());
1054              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1055                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1056              }
1057
1058              if ((*_blossom_data)[vb].status == MATCHED) {
1059                if (_blossom_set->classPrio(vb) ==
1060                    std::numeric_limits<Value>::max()) {
1061                  _delta2->erase(vb);
1062                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1063                           (*_blossom_data)[vb].offset) {
1064                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
1065                                   (*_blossom_data)[vb].offset);
1066                }
1067              }
1068            }
1069          }
1070        }
1071      }
1072    }
1073
1074    void oddToMatched(int blossom) {
1075      (*_blossom_data)[blossom].offset -= _delta_sum;
1076
1077      if (_blossom_set->classPrio(blossom) !=
1078          std::numeric_limits<Value>::max()) {
1079        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1080                       (*_blossom_data)[blossom].offset);
1081      }
1082
1083      if (!_blossom_set->trivial(blossom)) {
1084        _delta4->erase(blossom);
1085      }
1086    }
1087
1088    void oddToEven(int blossom, int tree) {
1089      if (!_blossom_set->trivial(blossom)) {
1090        _delta4->erase(blossom);
1091        (*_blossom_data)[blossom].pot -=
1092          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1093      }
1094
1095      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1096           n != INVALID; ++n) {
1097        int ni = (*_node_index)[n];
1098
1099        _blossom_set->increase(n, std::numeric_limits<Value>::max());
1100
1101        (*_node_data)[ni].heap.clear();
1102        (*_node_data)[ni].heap_index.clear();
1103        (*_node_data)[ni].pot +=
1104          2 * _delta_sum - (*_blossom_data)[blossom].offset;
1105
1106        _delta1->push(n, (*_node_data)[ni].pot);
1107
1108        for (InArcIt e(_graph, n); e != INVALID; ++e) {
1109          Node v = _graph.source(e);
1110          int vb = _blossom_set->find(v);
1111          int vi = (*_node_index)[v];
1112
1113          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1114            dualScale * _weight[e];
1115
1116          if ((*_blossom_data)[vb].status == EVEN) {
1117            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1118              _delta3->push(e, rw / 2);
1119            }
1120          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1121            if (_delta3->state(e) != _delta3->IN_HEAP) {
1122              _delta3->push(e, rw);
1123            }
1124          } else {
1125
1126            typename std::map<int, Arc>::iterator it =
1127              (*_node_data)[vi].heap_index.find(tree);
1128
1129            if (it != (*_node_data)[vi].heap_index.end()) {
1130              if ((*_node_data)[vi].heap[it->second] > rw) {
1131                (*_node_data)[vi].heap.replace(it->second, e);
1132                (*_node_data)[vi].heap.decrease(e, rw);
1133                it->second = e;
1134              }
1135            } else {
1136              (*_node_data)[vi].heap.push(e, rw);
1137              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1138            }
1139
1140            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1141              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1142
1143              if ((*_blossom_data)[vb].status == MATCHED) {
1144                if (_delta2->state(vb) != _delta2->IN_HEAP) {
1145                  _delta2->push(vb, _blossom_set->classPrio(vb) -
1146                               (*_blossom_data)[vb].offset);
1147                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1148                           (*_blossom_data)[vb].offset) {
1149                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1150                                   (*_blossom_data)[vb].offset);
1151                }
1152              }
1153            }
1154          }
1155        }
1156      }
1157      (*_blossom_data)[blossom].offset = 0;
1158    }
1159
1160
1161    void matchedToUnmatched(int blossom) {
1162      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1163        _delta2->erase(blossom);
1164      }
1165
1166      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1167           n != INVALID; ++n) {
1168        int ni = (*_node_index)[n];
1169
1170        _blossom_set->increase(n, std::numeric_limits<Value>::max());
1171
1172        (*_node_data)[ni].heap.clear();
1173        (*_node_data)[ni].heap_index.clear();
1174
1175        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1176          Node v = _graph.target(e);
1177          int vb = _blossom_set->find(v);
1178          int vi = (*_node_index)[v];
1179
1180          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1181            dualScale * _weight[e];
1182
1183          if ((*_blossom_data)[vb].status == EVEN) {
1184            if (_delta3->state(e) != _delta3->IN_HEAP) {
1185              _delta3->push(e, rw);
1186            }
1187          }
1188        }
1189      }
1190    }
1191
1192    void unmatchedToMatched(int blossom) {
1193      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1194           n != INVALID; ++n) {
1195        int ni = (*_node_index)[n];
1196
1197        for (InArcIt e(_graph, n); e != INVALID; ++e) {
1198          Node v = _graph.source(e);
1199          int vb = _blossom_set->find(v);
1200          int vi = (*_node_index)[v];
1201
1202          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1203            dualScale * _weight[e];
1204
1205          if (vb == blossom) {
1206            if (_delta3->state(e) == _delta3->IN_HEAP) {
1207              _delta3->erase(e);
1208            }
1209          } else if ((*_blossom_data)[vb].status == EVEN) {
1210
1211            if (_delta3->state(e) == _delta3->IN_HEAP) {
1212              _delta3->erase(e);
1213            }
1214
1215            int vt = _tree_set->find(vb);
1216
1217            Arc r = _graph.oppositeArc(e);
1218
1219            typename std::map<int, Arc>::iterator it =
1220              (*_node_data)[ni].heap_index.find(vt);
1221
1222            if (it != (*_node_data)[ni].heap_index.end()) {
1223              if ((*_node_data)[ni].heap[it->second] > rw) {
1224                (*_node_data)[ni].heap.replace(it->second, r);
1225                (*_node_data)[ni].heap.decrease(r, rw);
1226                it->second = r;
1227              }
1228            } else {
1229              (*_node_data)[ni].heap.push(r, rw);
1230              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1231            }
1232
1233            if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1234              _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1235
1236              if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1237                _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1238                             (*_blossom_data)[blossom].offset);
1239              } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
1240                         (*_blossom_data)[blossom].offset){
1241                _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1242                                 (*_blossom_data)[blossom].offset);
1243              }
1244            }
1245
1246          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1247            if (_delta3->state(e) == _delta3->IN_HEAP) {
1248              _delta3->erase(e);
1249            }
1250          }
1251        }
1252      }
1253    }
1254
1255    void alternatePath(int even, int tree) {
1256      int odd;
1257
1258      evenToMatched(even, tree);
1259      (*_blossom_data)[even].status = MATCHED;
1260
1261      while ((*_blossom_data)[even].pred != INVALID) {
1262        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
1263        (*_blossom_data)[odd].status = MATCHED;
1264        oddToMatched(odd);
1265        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1266
1267        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
1268        (*_blossom_data)[even].status = MATCHED;
1269        evenToMatched(even, tree);
1270        (*_blossom_data)[even].next =
1271          _graph.oppositeArc((*_blossom_data)[odd].pred);
1272      }
1273
1274    }
1275
1276    void destroyTree(int tree) {
1277      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1278        if ((*_blossom_data)[b].status == EVEN) {
1279          (*_blossom_data)[b].status = MATCHED;
1280          evenToMatched(b, tree);
1281        } else if ((*_blossom_data)[b].status == ODD) {
1282          (*_blossom_data)[b].status = MATCHED;
1283          oddToMatched(b);
1284        }
1285      }
1286      _tree_set->eraseClass(tree);
1287    }
1288
1289
1290    void unmatchNode(const Node& node) {
1291      int blossom = _blossom_set->find(node);
1292      int tree = _tree_set->find(blossom);
1293
1294      alternatePath(blossom, tree);
1295      destroyTree(tree);
1296
1297      (*_blossom_data)[blossom].status = UNMATCHED;
1298      (*_blossom_data)[blossom].base = node;
1299      matchedToUnmatched(blossom);
1300    }
1301
1302
1303    void augmentOnEdge(const Edge& edge) {
1304
1305      int left = _blossom_set->find(_graph.u(edge));
1306      int right = _blossom_set->find(_graph.v(edge));
1307
1308      if ((*_blossom_data)[left].status == EVEN) {
1309        int left_tree = _tree_set->find(left);
1310        alternatePath(left, left_tree);
1311        destroyTree(left_tree);
1312      } else {
1313        (*_blossom_data)[left].status = MATCHED;
1314        unmatchedToMatched(left);
1315      }
1316
1317      if ((*_blossom_data)[right].status == EVEN) {
1318        int right_tree = _tree_set->find(right);
1319        alternatePath(right, right_tree);
1320        destroyTree(right_tree);
1321      } else {
1322        (*_blossom_data)[right].status = MATCHED;
1323        unmatchedToMatched(right);
1324      }
1325
1326      (*_blossom_data)[left].next = _graph.direct(edge, true);
1327      (*_blossom_data)[right].next = _graph.direct(edge, false);
1328    }
1329
1330    void extendOnArc(const Arc& arc) {
1331      int base = _blossom_set->find(_graph.target(arc));
1332      int tree = _tree_set->find(base);
1333
1334      int odd = _blossom_set->find(_graph.source(arc));
1335      _tree_set->insert(odd, tree);
1336      (*_blossom_data)[odd].status = ODD;
1337      matchedToOdd(odd);
1338      (*_blossom_data)[odd].pred = arc;
1339
1340      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
1341      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1342      _tree_set->insert(even, tree);
1343      (*_blossom_data)[even].status = EVEN;
1344      matchedToEven(even, tree);
1345    }
1346
1347    void shrinkOnEdge(const Edge& edge, int tree) {
1348      int nca = -1;
1349      std::vector<int> left_path, right_path;
1350
1351      {
1352        std::set<int> left_set, right_set;
1353        int left = _blossom_set->find(_graph.u(edge));
1354        left_path.push_back(left);
1355        left_set.insert(left);
1356
1357        int right = _blossom_set->find(_graph.v(edge));
1358        right_path.push_back(right);
1359        right_set.insert(right);
1360
1361        while (true) {
1362
1363          if ((*_blossom_data)[left].pred == INVALID) break;
1364
1365          left =
1366            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1367          left_path.push_back(left);
1368          left =
1369            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1370          left_path.push_back(left);
1371
1372          left_set.insert(left);
1373
1374          if (right_set.find(left) != right_set.end()) {
1375            nca = left;
1376            break;
1377          }
1378
1379          if ((*_blossom_data)[right].pred == INVALID) break;
1380
1381          right =
1382            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1383          right_path.push_back(right);
1384          right =
1385            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1386          right_path.push_back(right);
1387
1388          right_set.insert(right);
1389
1390          if (left_set.find(right) != left_set.end()) {
1391            nca = right;
1392            break;
1393          }
1394
1395        }
1396
1397        if (nca == -1) {
1398          if ((*_blossom_data)[left].pred == INVALID) {
1399            nca = right;
1400            while (left_set.find(nca) == left_set.end()) {
1401              nca =
1402                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1403              right_path.push_back(nca);
1404              nca =
1405                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1406              right_path.push_back(nca);
1407            }
1408          } else {
1409            nca = left;
1410            while (right_set.find(nca) == right_set.end()) {
1411              nca =
1412                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1413              left_path.push_back(nca);
1414              nca =
1415                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1416              left_path.push_back(nca);
1417            }
1418          }
1419        }
1420      }
1421
1422      std::vector<int> subblossoms;
1423      Arc prev;
1424
1425      prev = _graph.direct(edge, true);
1426      for (int i = 0; left_path[i] != nca; i += 2) {
1427        subblossoms.push_back(left_path[i]);
1428        (*_blossom_data)[left_path[i]].next = prev;
1429        _tree_set->erase(left_path[i]);
1430
1431        subblossoms.push_back(left_path[i + 1]);
1432        (*_blossom_data)[left_path[i + 1]].status = EVEN;
1433        oddToEven(left_path[i + 1], tree);
1434        _tree_set->erase(left_path[i + 1]);
1435        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
1436      }
1437
1438      int k = 0;
1439      while (right_path[k] != nca) ++k;
1440
1441      subblossoms.push_back(nca);
1442      (*_blossom_data)[nca].next = prev;
1443
1444      for (int i = k - 2; i >= 0; i -= 2) {
1445        subblossoms.push_back(right_path[i + 1]);
1446        (*_blossom_data)[right_path[i + 1]].status = EVEN;
1447        oddToEven(right_path[i + 1], tree);
1448        _tree_set->erase(right_path[i + 1]);
1449
1450        (*_blossom_data)[right_path[i + 1]].next =
1451          (*_blossom_data)[right_path[i + 1]].pred;
1452
1453        subblossoms.push_back(right_path[i]);
1454        _tree_set->erase(right_path[i]);
1455      }
1456
1457      int surface =
1458        _blossom_set->join(subblossoms.begin(), subblossoms.end());
1459
1460      for (int i = 0; i < int(subblossoms.size()); ++i) {
1461        if (!_blossom_set->trivial(subblossoms[i])) {
1462          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1463        }
1464        (*_blossom_data)[subblossoms[i]].status = MATCHED;
1465      }
1466
1467      (*_blossom_data)[surface].pot = -2 * _delta_sum;
1468      (*_blossom_data)[surface].offset = 0;
1469      (*_blossom_data)[surface].status = EVEN;
1470      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1471      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1472
1473      _tree_set->insert(surface, tree);
1474      _tree_set->erase(nca);
1475    }
1476
1477    void splitBlossom(int blossom) {
1478      Arc next = (*_blossom_data)[blossom].next;
1479      Arc pred = (*_blossom_data)[blossom].pred;
1480
1481      int tree = _tree_set->find(blossom);
1482
1483      (*_blossom_data)[blossom].status = MATCHED;
1484      oddToMatched(blossom);
1485      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1486        _delta2->erase(blossom);
1487      }
1488
1489      std::vector<int> subblossoms;
1490      _blossom_set->split(blossom, std::back_inserter(subblossoms));
1491
1492      Value offset = (*_blossom_data)[blossom].offset;
1493      int b = _blossom_set->find(_graph.source(pred));
1494      int d = _blossom_set->find(_graph.source(next));
1495
1496      int ib = -1, id = -1;
1497      for (int i = 0; i < int(subblossoms.size()); ++i) {
1498        if (subblossoms[i] == b) ib = i;
1499        if (subblossoms[i] == d) id = i;
1500
1501        (*_blossom_data)[subblossoms[i]].offset = offset;
1502        if (!_blossom_set->trivial(subblossoms[i])) {
1503          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1504        }
1505        if (_blossom_set->classPrio(subblossoms[i]) !=
1506            std::numeric_limits<Value>::max()) {
1507          _delta2->push(subblossoms[i],
1508                        _blossom_set->classPrio(subblossoms[i]) -
1509                        (*_blossom_data)[subblossoms[i]].offset);
1510        }
1511      }
1512
1513      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1514        for (int i = (id + 1) % subblossoms.size();
1515             i != ib; i = (i + 2) % subblossoms.size()) {
1516          int sb = subblossoms[i];
1517          int tb = subblossoms[(i + 1) % subblossoms.size()];
1518          (*_blossom_data)[sb].next =
1519            _graph.oppositeArc((*_blossom_data)[tb].next);
1520        }
1521
1522        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1523          int sb = subblossoms[i];
1524          int tb = subblossoms[(i + 1) % subblossoms.size()];
1525          int ub = subblossoms[(i + 2) % subblossoms.size()];
1526
1527          (*_blossom_data)[sb].status = ODD;
1528          matchedToOdd(sb);
1529          _tree_set->insert(sb, tree);
1530          (*_blossom_data)[sb].pred = pred;
1531          (*_blossom_data)[sb].next =
1532                           _graph.oppositeArc((*_blossom_data)[tb].next);
1533
1534          pred = (*_blossom_data)[ub].next;
1535
1536          (*_blossom_data)[tb].status = EVEN;
1537          matchedToEven(tb, tree);
1538          _tree_set->insert(tb, tree);
1539          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1540        }
1541
1542        (*_blossom_data)[subblossoms[id]].status = ODD;
1543        matchedToOdd(subblossoms[id]);
1544        _tree_set->insert(subblossoms[id], tree);
1545        (*_blossom_data)[subblossoms[id]].next = next;
1546        (*_blossom_data)[subblossoms[id]].pred = pred;
1547
1548      } else {
1549
1550        for (int i = (ib + 1) % subblossoms.size();
1551             i != id; i = (i + 2) % subblossoms.size()) {
1552          int sb = subblossoms[i];
1553          int tb = subblossoms[(i + 1) % subblossoms.size()];
1554          (*_blossom_data)[sb].next =
1555            _graph.oppositeArc((*_blossom_data)[tb].next);
1556        }
1557
1558        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1559          int sb = subblossoms[i];
1560          int tb = subblossoms[(i + 1) % subblossoms.size()];
1561          int ub = subblossoms[(i + 2) % subblossoms.size()];
1562
1563          (*_blossom_data)[sb].status = ODD;
1564          matchedToOdd(sb);
1565          _tree_set->insert(sb, tree);
1566          (*_blossom_data)[sb].next = next;
1567          (*_blossom_data)[sb].pred =
1568            _graph.oppositeArc((*_blossom_data)[tb].next);
1569
1570          (*_blossom_data)[tb].status = EVEN;
1571          matchedToEven(tb, tree);
1572          _tree_set->insert(tb, tree);
1573          (*_blossom_data)[tb].pred =
1574            (*_blossom_data)[tb].next =
1575            _graph.oppositeArc((*_blossom_data)[ub].next);
1576          next = (*_blossom_data)[ub].next;
1577        }
1578
1579        (*_blossom_data)[subblossoms[ib]].status = ODD;
1580        matchedToOdd(subblossoms[ib]);
1581        _tree_set->insert(subblossoms[ib], tree);
1582        (*_blossom_data)[subblossoms[ib]].next = next;
1583        (*_blossom_data)[subblossoms[ib]].pred = pred;
1584      }
1585      _tree_set->erase(blossom);
1586    }
1587
1588    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1589      if (_blossom_set->trivial(blossom)) {
1590        int bi = (*_node_index)[base];
1591        Value pot = (*_node_data)[bi].pot;
1592
1593        (*_matching)[base] = matching;
1594        _blossom_node_list.push_back(base);
1595        (*_node_potential)[base] = pot;
1596      } else {
1597
1598        Value pot = (*_blossom_data)[blossom].pot;
1599        int bn = _blossom_node_list.size();
1600
1601        std::vector<int> subblossoms;
1602        _blossom_set->split(blossom, std::back_inserter(subblossoms));
1603        int b = _blossom_set->find(base);
1604        int ib = -1;
1605        for (int i = 0; i < int(subblossoms.size()); ++i) {
1606          if (subblossoms[i] == b) { ib = i; break; }
1607        }
1608
1609        for (int i = 1; i < int(subblossoms.size()); i += 2) {
1610          int sb = subblossoms[(ib + i) % subblossoms.size()];
1611          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1612
1613          Arc m = (*_blossom_data)[tb].next;
1614          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1615          extractBlossom(tb, _graph.source(m), m);
1616        }
1617        extractBlossom(subblossoms[ib], base, matching);
1618
1619        int en = _blossom_node_list.size();
1620
1621        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1622      }
1623    }
1624
1625    void extractMatching() {
1626      std::vector<int> blossoms;
1627      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1628        blossoms.push_back(c);
1629      }
1630
1631      for (int i = 0; i < int(blossoms.size()); ++i) {
1632        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
1633
1634          Value offset = (*_blossom_data)[blossoms[i]].offset;
1635          (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1636          for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1637               n != INVALID; ++n) {
1638            (*_node_data)[(*_node_index)[n]].pot -= offset;
1639          }
1640
1641          Arc matching = (*_blossom_data)[blossoms[i]].next;
1642          Node base = _graph.source(matching);
1643          extractBlossom(blossoms[i], base, matching);
1644        } else {
1645          Node base = (*_blossom_data)[blossoms[i]].base;
1646          extractBlossom(blossoms[i], base, INVALID);
1647        }
1648      }
1649    }
1650
1651  public:
1652
1653    /// \brief Constructor
1654    ///
1655    /// Constructor.
1656    MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1657      : _graph(graph), _weight(weight), _matching(0),
1658        _node_potential(0), _blossom_potential(), _blossom_node_list(),
1659        _node_num(0), _blossom_num(0),
1660
1661        _blossom_index(0), _blossom_set(0), _blossom_data(0),
1662        _node_index(0), _node_heap_index(0), _node_data(0),
1663        _tree_set_index(0), _tree_set(0),
1664
1665        _delta1_index(0), _delta1(0),
1666        _delta2_index(0), _delta2(0),
1667        _delta3_index(0), _delta3(0),
1668        _delta4_index(0), _delta4(0),
1669
1670        _delta_sum() {}
1671
1672    ~MaxWeightedMatching() {
1673      destroyStructures();
1674    }
1675
1676    /// \name Execution Control
1677    /// The simplest way to execute the algorithm is to use the
1678    /// \ref run() member function.
1679
1680    ///@{
1681
1682    /// \brief Initialize the algorithm
1683    ///
1684    /// This function initializes the algorithm.
1685    void init() {
1686      createStructures();
1687
1688      for (ArcIt e(_graph); e != INVALID; ++e) {
1689        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1690      }
1691      for (NodeIt n(_graph); n != INVALID; ++n) {
1692        (*_delta1_index)[n] = _delta1->PRE_HEAP;
1693      }
1694      for (EdgeIt e(_graph); e != INVALID; ++e) {
1695        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1696      }
1697      for (int i = 0; i < _blossom_num; ++i) {
1698        (*_delta2_index)[i] = _delta2->PRE_HEAP;
1699        (*_delta4_index)[i] = _delta4->PRE_HEAP;
1700      }
1701
1702      int index = 0;
1703      for (NodeIt n(_graph); n != INVALID; ++n) {
1704        Value max = 0;
1705        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1706          if (_graph.target(e) == n) continue;
1707          if ((dualScale * _weight[e]) / 2 > max) {
1708            max = (dualScale * _weight[e]) / 2;
1709          }
1710        }
1711        (*_node_index)[n] = index;
1712        (*_node_data)[index].pot = max;
1713        _delta1->push(n, max);
1714        int blossom =
1715          _blossom_set->insert(n, std::numeric_limits<Value>::max());
1716
1717        _tree_set->insert(blossom);
1718
1719        (*_blossom_data)[blossom].status = EVEN;
1720        (*_blossom_data)[blossom].pred = INVALID;
1721        (*_blossom_data)[blossom].next = INVALID;
1722        (*_blossom_data)[blossom].pot = 0;
1723        (*_blossom_data)[blossom].offset = 0;
1724        ++index;
1725      }
1726      for (EdgeIt e(_graph); e != INVALID; ++e) {
1727        int si = (*_node_index)[_graph.u(e)];
1728        int ti = (*_node_index)[_graph.v(e)];
1729        if (_graph.u(e) != _graph.v(e)) {
1730          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1731                            dualScale * _weight[e]) / 2);
1732        }
1733      }
1734    }
1735
1736    /// \brief Start the algorithm
1737    ///
1738    /// This function starts the algorithm.
1739    ///
1740    /// \pre \ref init() must be called before using this function.
1741    void start() {
1742      enum OpType {
1743        D1, D2, D3, D4
1744      };
1745
1746      int unmatched = _node_num;
1747      while (unmatched > 0) {
1748        Value d1 = !_delta1->empty() ?
1749          _delta1->prio() : std::numeric_limits<Value>::max();
1750
1751        Value d2 = !_delta2->empty() ?
1752          _delta2->prio() : std::numeric_limits<Value>::max();
1753
1754        Value d3 = !_delta3->empty() ?
1755          _delta3->prio() : std::numeric_limits<Value>::max();
1756
1757        Value d4 = !_delta4->empty() ?
1758          _delta4->prio() : std::numeric_limits<Value>::max();
1759
1760        _delta_sum = d1; OpType ot = D1;
1761        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1762        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1763        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1764
1765
1766        switch (ot) {
1767        case D1:
1768          {
1769            Node n = _delta1->top();
1770            unmatchNode(n);
1771            --unmatched;
1772          }
1773          break;
1774        case D2:
1775          {
1776            int blossom = _delta2->top();
1777            Node n = _blossom_set->classTop(blossom);
1778            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
1779            extendOnArc(e);
1780          }
1781          break;
1782        case D3:
1783          {
1784            Edge e = _delta3->top();
1785
1786            int left_blossom = _blossom_set->find(_graph.u(e));
1787            int right_blossom = _blossom_set->find(_graph.v(e));
1788
1789            if (left_blossom == right_blossom) {
1790              _delta3->pop();
1791            } else {
1792              int left_tree;
1793              if ((*_blossom_data)[left_blossom].status == EVEN) {
1794                left_tree = _tree_set->find(left_blossom);
1795              } else {
1796                left_tree = -1;
1797                ++unmatched;
1798              }
1799              int right_tree;
1800              if ((*_blossom_data)[right_blossom].status == EVEN) {
1801                right_tree = _tree_set->find(right_blossom);
1802              } else {
1803                right_tree = -1;
1804                ++unmatched;
1805              }
1806
1807              if (left_tree == right_tree) {
1808                shrinkOnEdge(e, left_tree);
1809              } else {
1810                augmentOnEdge(e);
1811                unmatched -= 2;
1812              }
1813            }
1814          } break;
1815        case D4:
1816          splitBlossom(_delta4->top());
1817          break;
1818        }
1819      }
1820      extractMatching();
1821    }
1822
1823    /// \brief Run the algorithm.
1824    ///
1825    /// This method runs the \c %MaxWeightedMatching algorithm.
1826    ///
1827    /// \note mwm.run() is just a shortcut of the following code.
1828    /// \code
1829    ///   mwm.init();
1830    ///   mwm.start();
1831    /// \endcode
1832    void run() {
1833      init();
1834      start();
1835    }
1836
1837    /// @}
1838
1839    /// \name Primal Solution
1840    /// Functions to get the primal solution, i.e. the maximum weighted
1841    /// matching.\n
1842    /// Either \ref run() or \ref start() function should be called before
1843    /// using them.
1844
1845    /// @{
1846
1847    /// \brief Return the weight of the matching.
1848    ///
1849    /// This function returns the weight of the found matching.
1850    ///
1851    /// \pre Either run() or start() must be called before using this function.
1852    Value matchingWeight() const {
1853      Value sum = 0;
1854      for (NodeIt n(_graph); n != INVALID; ++n) {
1855        if ((*_matching)[n] != INVALID) {
1856          sum += _weight[(*_matching)[n]];
1857        }
1858      }
1859      return sum /= 2;
1860    }
1861
1862    /// \brief Return the size (cardinality) of the matching.
1863    ///
1864    /// This function returns the size (cardinality) of the found matching.
1865    ///
1866    /// \pre Either run() or start() must be called before using this function.
1867    int matchingSize() const {
1868      int num = 0;
1869      for (NodeIt n(_graph); n != INVALID; ++n) {
1870        if ((*_matching)[n] != INVALID) {
1871          ++num;
1872        }
1873      }
1874      return num /= 2;
1875    }
1876
1877    /// \brief Return \c true if the given edge is in the matching.
1878    ///
1879    /// This function returns \c true if the given edge is in the found
1880    /// matching.
1881    ///
1882    /// \pre Either run() or start() must be called before using this function.
1883    bool matching(const Edge& edge) const {
1884      return edge == (*_matching)[_graph.u(edge)];
1885    }
1886
1887    /// \brief Return the matching arc (or edge) incident to the given node.
1888    ///
1889    /// This function returns the matching arc (or edge) incident to the
1890    /// given node in the found matching or \c INVALID if the node is
1891    /// not covered by the matching.
1892    ///
1893    /// \pre Either run() or start() must be called before using this function.
1894    Arc matching(const Node& node) const {
1895      return (*_matching)[node];
1896    }
1897
1898    /// \brief Return a const reference to the matching map.
1899    ///
1900    /// This function returns a const reference to a node map that stores
1901    /// the matching arc (or edge) incident to each node.
1902    const MatchingMap& matchingMap() const {
1903      return *_matching;
1904    }
1905
1906    /// \brief Return the mate of the given node.
1907    ///
1908    /// This function returns the mate of the given node in the found
1909    /// matching or \c INVALID if the node is not covered by the matching.
1910    ///
1911    /// \pre Either run() or start() must be called before using this function.
1912    Node mate(const Node& node) const {
1913      return (*_matching)[node] != INVALID ?
1914        _graph.target((*_matching)[node]) : INVALID;
1915    }
1916
1917    /// @}
1918
1919    /// \name Dual Solution
1920    /// Functions to get the dual solution.\n
1921    /// Either \ref run() or \ref start() function should be called before
1922    /// using them.
1923
1924    /// @{
1925
1926    /// \brief Return the value of the dual solution.
1927    ///
1928    /// This function returns the value of the dual solution.
1929    /// It should be equal to the primal value scaled by \ref dualScale
1930    /// "dual scale".
1931    ///
1932    /// \pre Either run() or start() must be called before using this function.
1933    Value dualValue() const {
1934      Value sum = 0;
1935      for (NodeIt n(_graph); n != INVALID; ++n) {
1936        sum += nodeValue(n);
1937      }
1938      for (int i = 0; i < blossomNum(); ++i) {
1939        sum += blossomValue(i) * (blossomSize(i) / 2);
1940      }
1941      return sum;
1942    }
1943
1944    /// \brief Return the dual value (potential) of the given node.
1945    ///
1946    /// This function returns the dual value (potential) of the given node.
1947    ///
1948    /// \pre Either run() or start() must be called before using this function.
1949    Value nodeValue(const Node& n) const {
1950      return (*_node_potential)[n];
1951    }
1952
1953    /// \brief Return the number of the blossoms in the basis.
1954    ///
1955    /// This function returns the number of the blossoms in the basis.
1956    ///
1957    /// \pre Either run() or start() must be called before using this function.
1958    /// \see BlossomIt
1959    int blossomNum() const {
1960      return _blossom_potential.size();
1961    }
1962
1963    /// \brief Return the number of the nodes in the given blossom.
1964    ///
1965    /// This function returns the number of the nodes in the given blossom.
1966    ///
1967    /// \pre Either run() or start() must be called before using this function.
1968    /// \see BlossomIt
1969    int blossomSize(int k) const {
1970      return _blossom_potential[k].end - _blossom_potential[k].begin;
1971    }
1972
1973    /// \brief Return the dual value (ptential) of the given blossom.
1974    ///
1975    /// This function returns the dual value (ptential) of the given blossom.
1976    ///
1977    /// \pre Either run() or start() must be called before using this function.
1978    Value blossomValue(int k) const {
1979      return _blossom_potential[k].value;
1980    }
1981
1982    /// \brief Iterator for obtaining the nodes of a blossom.
1983    ///
1984    /// This class provides an iterator for obtaining the nodes of the
1985    /// given blossom. It lists a subset of the nodes.
1986    /// Before using this iterator, you must allocate a
1987    /// MaxWeightedMatching class and execute it.
1988    class BlossomIt {
1989    public:
1990
1991      /// \brief Constructor.
1992      ///
1993      /// Constructor to get the nodes of the given variable.
1994      ///
1995      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
1996      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
1997      /// called before initializing this iterator.
1998      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
1999        : _algorithm(&algorithm)
2000      {
2001        _index = _algorithm->_blossom_potential[variable].begin;
2002        _last = _algorithm->_blossom_potential[variable].end;
2003      }
2004
2005      /// \brief Conversion to \c Node.
2006      ///
2007      /// Conversion to \c Node.
2008      operator Node() const {
2009        return _algorithm->_blossom_node_list[_index];
2010      }
2011
2012      /// \brief Increment operator.
2013      ///
2014      /// Increment operator.
2015      BlossomIt& operator++() {
2016        ++_index;
2017        return *this;
2018      }
2019
2020      /// \brief Validity checking
2021      ///
2022      /// Checks whether the iterator is invalid.
2023      bool operator==(Invalid) const { return _index == _last; }
2024
2025      /// \brief Validity checking
2026      ///
2027      /// Checks whether the iterator is valid.
2028      bool operator!=(Invalid) const { return _index != _last; }
2029
2030    private:
2031      const MaxWeightedMatching* _algorithm;
2032      int _last;
2033      int _index;
2034    };
2035
2036    /// @}
2037
2038  };
2039
2040  /// \ingroup matching
2041  ///
2042  /// \brief Weighted perfect matching in general graphs
2043  ///
2044  /// This class provides an efficient implementation of Edmond's
2045  /// maximum weighted perfect matching algorithm. The implementation
2046  /// is based on extensive use of priority queues and provides
2047  /// \f$O(nm\log n)\f$ time complexity.
2048  ///
2049  /// The maximum weighted perfect matching problem is to find a subset of
2050  /// the edges in an undirected graph with maximum overall weight for which
2051  /// each node has exactly one incident edge.
2052  /// It can be formulated with the following linear program.
2053  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
2054  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
2055      \quad \forall B\in\mathcal{O}\f] */
2056  /// \f[x_e \ge 0\quad \forall e\in E\f]
2057  /// \f[\max \sum_{e\in E}x_ew_e\f]
2058  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
2059  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
2060  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
2061  /// subsets of the nodes.
2062  ///
2063  /// The algorithm calculates an optimal matching and a proof of the
2064  /// optimality. The solution of the dual problem can be used to check
2065  /// the result of the algorithm. The dual linear problem is the
2066  /// following.
2067  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
2068      w_{uv} \quad \forall uv\in E\f] */
2069  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
2070  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
2071      \frac{\vert B \vert - 1}{2}z_B\f] */
2072  ///
2073  /// The algorithm can be executed with the run() function.
2074  /// After it the matching (the primal solution) and the dual solution
2075  /// can be obtained using the query functions and the
2076  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
2077  /// which is able to iterate on the nodes of a blossom.
2078  /// If the value type is integer, then the dual solution is multiplied
2079  /// by \ref MaxWeightedMatching::dualScale "4".
2080  ///
2081  /// \tparam GR The undirected graph type the algorithm runs on.
2082  /// \tparam WM The type edge weight map. The default type is
2083  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
2084#ifdef DOXYGEN
2085  template <typename GR, typename WM>
2086#else
2087  template <typename GR,
2088            typename WM = typename GR::template EdgeMap<int> >
2089#endif
2090  class MaxWeightedPerfectMatching {
2091  public:
2092
2093    /// The graph type of the algorithm
2094    typedef GR Graph;
2095    /// The type of the edge weight map
2096    typedef WM WeightMap;
2097    /// The value type of the edge weights
2098    typedef typename WeightMap::Value Value;
2099
2100    /// \brief Scaling factor for dual solution
2101    ///
2102    /// Scaling factor for dual solution, it is equal to 4 or 1
2103    /// according to the value type.
2104    static const int dualScale =
2105      std::numeric_limits<Value>::is_integer ? 4 : 1;
2106
2107    /// The type of the matching map
2108    typedef typename Graph::template NodeMap<typename Graph::Arc>
2109    MatchingMap;
2110
2111  private:
2112
2113    TEMPLATE_GRAPH_TYPEDEFS(Graph);
2114
2115    typedef typename Graph::template NodeMap<Value> NodePotential;
2116    typedef std::vector<Node> BlossomNodeList;
2117
2118    struct BlossomVariable {
2119      int begin, end;
2120      Value value;
2121
2122      BlossomVariable(int _begin, int _end, Value _value)
2123        : begin(_begin), end(_end), value(_value) {}
2124
2125    };
2126
2127    typedef std::vector<BlossomVariable> BlossomPotential;
2128
2129    const Graph& _graph;
2130    const WeightMap& _weight;
2131
2132    MatchingMap* _matching;
2133
2134    NodePotential* _node_potential;
2135
2136    BlossomPotential _blossom_potential;
2137    BlossomNodeList _blossom_node_list;
2138
2139    int _node_num;
2140    int _blossom_num;
2141
2142    typedef RangeMap<int> IntIntMap;
2143
2144    enum Status {
2145      EVEN = -1, MATCHED = 0, ODD = 1
2146    };
2147
2148    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
2149    struct BlossomData {
2150      int tree;
2151      Status status;
2152      Arc pred, next;
2153      Value pot, offset;
2154    };
2155
2156    IntNodeMap *_blossom_index;
2157    BlossomSet *_blossom_set;
2158    RangeMap<BlossomData>* _blossom_data;
2159
2160    IntNodeMap *_node_index;
2161    IntArcMap *_node_heap_index;
2162
2163    struct NodeData {
2164
2165      NodeData(IntArcMap& node_heap_index)
2166        : heap(node_heap_index) {}
2167
2168      int blossom;
2169      Value pot;
2170      BinHeap<Value, IntArcMap> heap;
2171      std::map<int, Arc> heap_index;
2172
2173      int tree;
2174    };
2175
2176    RangeMap<NodeData>* _node_data;
2177
2178    typedef ExtendFindEnum<IntIntMap> TreeSet;
2179
2180    IntIntMap *_tree_set_index;
2181    TreeSet *_tree_set;
2182
2183    IntIntMap *_delta2_index;
2184    BinHeap<Value, IntIntMap> *_delta2;
2185
2186    IntEdgeMap *_delta3_index;
2187    BinHeap<Value, IntEdgeMap> *_delta3;
2188
2189    IntIntMap *_delta4_index;
2190    BinHeap<Value, IntIntMap> *_delta4;
2191
2192    Value _delta_sum;
2193
2194    void createStructures() {
2195      _node_num = countNodes(_graph);
2196      _blossom_num = _node_num * 3 / 2;
2197
2198      if (!_matching) {
2199        _matching = new MatchingMap(_graph);
2200      }
2201      if (!_node_potential) {
2202        _node_potential = new NodePotential(_graph);
2203      }
2204      if (!_blossom_set) {
2205        _blossom_index = new IntNodeMap(_graph);
2206        _blossom_set = new BlossomSet(*_blossom_index);
2207        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2208      }
2209
2210      if (!_node_index) {
2211        _node_index = new IntNodeMap(_graph);
2212        _node_heap_index = new IntArcMap(_graph);
2213        _node_data = new RangeMap<NodeData>(_node_num,
2214                                            NodeData(*_node_heap_index));
2215      }
2216
2217      if (!_tree_set) {
2218        _tree_set_index = new IntIntMap(_blossom_num);
2219        _tree_set = new TreeSet(*_tree_set_index);
2220      }
2221      if (!_delta2) {
2222        _delta2_index = new IntIntMap(_blossom_num);
2223        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2224      }
2225      if (!_delta3) {
2226        _delta3_index = new IntEdgeMap(_graph);
2227        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
2228      }
2229      if (!_delta4) {
2230        _delta4_index = new IntIntMap(_blossom_num);
2231        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2232      }
2233    }
2234
2235    void destroyStructures() {
2236      _node_num = countNodes(_graph);
2237      _blossom_num = _node_num * 3 / 2;
2238
2239      if (_matching) {
2240        delete _matching;
2241      }
2242      if (_node_potential) {
2243        delete _node_potential;
2244      }
2245      if (_blossom_set) {
2246        delete _blossom_index;
2247        delete _blossom_set;
2248        delete _blossom_data;
2249      }
2250
2251      if (_node_index) {
2252        delete _node_index;
2253        delete _node_heap_index;
2254        delete _node_data;
2255      }
2256
2257      if (_tree_set) {
2258        delete _tree_set_index;
2259        delete _tree_set;
2260      }
2261      if (_delta2) {
2262        delete _delta2_index;
2263        delete _delta2;
2264      }
2265      if (_delta3) {
2266        delete _delta3_index;
2267        delete _delta3;
2268      }
2269      if (_delta4) {
2270        delete _delta4_index;
2271        delete _delta4;
2272      }
2273    }
2274
2275    void matchedToEven(int blossom, int tree) {
2276      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2277        _delta2->erase(blossom);
2278      }
2279
2280      if (!_blossom_set->trivial(blossom)) {
2281        (*_blossom_data)[blossom].pot -=
2282          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2283      }
2284
2285      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2286           n != INVALID; ++n) {
2287
2288        _blossom_set->increase(n, std::numeric_limits<Value>::max());
2289        int ni = (*_node_index)[n];
2290
2291        (*_node_data)[ni].heap.clear();
2292        (*_node_data)[ni].heap_index.clear();
2293
2294        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
2295
2296        for (InArcIt e(_graph, n); e != INVALID; ++e) {
2297          Node v = _graph.source(e);
2298          int vb = _blossom_set->find(v);
2299          int vi = (*_node_index)[v];
2300
2301          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2302            dualScale * _weight[e];
2303
2304          if ((*_blossom_data)[vb].status == EVEN) {
2305            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2306              _delta3->push(e, rw / 2);
2307            }
2308          } else {
2309            typename std::map<int, Arc>::iterator it =
2310              (*_node_data)[vi].heap_index.find(tree);
2311
2312            if (it != (*_node_data)[vi].heap_index.end()) {
2313              if ((*_node_data)[vi].heap[it->second] > rw) {
2314                (*_node_data)[vi].heap.replace(it->second, e);
2315                (*_node_data)[vi].heap.decrease(e, rw);
2316                it->second = e;
2317              }
2318            } else {
2319              (*_node_data)[vi].heap.push(e, rw);
2320              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2321            }
2322
2323            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2324              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2325
2326              if ((*_blossom_data)[vb].status == MATCHED) {
2327                if (_delta2->state(vb) != _delta2->IN_HEAP) {
2328                  _delta2->push(vb, _blossom_set->classPrio(vb) -
2329                               (*_blossom_data)[vb].offset);
2330                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2331                           (*_blossom_data)[vb].offset){
2332                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2333                                   (*_blossom_data)[vb].offset);
2334                }
2335              }
2336            }
2337          }
2338        }
2339      }
2340      (*_blossom_data)[blossom].offset = 0;
2341    }
2342
2343    void matchedToOdd(int blossom) {
2344      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2345        _delta2->erase(blossom);
2346      }
2347      (*_blossom_data)[blossom].offset += _delta_sum;
2348      if (!_blossom_set->trivial(blossom)) {
2349        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
2350                     (*_blossom_data)[blossom].offset);
2351      }
2352    }
2353
2354    void evenToMatched(int blossom, int tree) {
2355      if (!_blossom_set->trivial(blossom)) {
2356        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
2357      }
2358
2359      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2360           n != INVALID; ++n) {
2361        int ni = (*_node_index)[n];
2362        (*_node_data)[ni].pot -= _delta_sum;
2363
2364        for (InArcIt e(_graph, n); e != INVALID; ++e) {
2365          Node v = _graph.source(e);
2366          int vb = _blossom_set->find(v);
2367          int vi = (*_node_index)[v];
2368
2369          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2370            dualScale * _weight[e];
2371
2372          if (vb == blossom) {
2373            if (_delta3->state(e) == _delta3->IN_HEAP) {
2374              _delta3->erase(e);
2375            }
2376          } else if ((*_blossom_data)[vb].status == EVEN) {
2377
2378            if (_delta3->state(e) == _delta3->IN_HEAP) {
2379              _delta3->erase(e);
2380            }
2381
2382            int vt = _tree_set->find(vb);
2383
2384            if (vt != tree) {
2385
2386              Arc r = _graph.oppositeArc(e);
2387
2388              typename std::map<int, Arc>::iterator it =
2389                (*_node_data)[ni].heap_index.find(vt);
2390
2391              if (it != (*_node_data)[ni].heap_index.end()) {
2392                if ((*_node_data)[ni].heap[it->second] > rw) {
2393                  (*_node_data)[ni].heap.replace(it->second, r);
2394                  (*_node_data)[ni].heap.decrease(r, rw);
2395                  it->second = r;
2396                }
2397              } else {
2398                (*_node_data)[ni].heap.push(r, rw);
2399                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2400              }
2401
2402              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2403                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2404
2405                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2406                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2407                               (*_blossom_data)[blossom].offset);
2408                } else if ((*_delta2)[blossom] >
2409                           _blossom_set->classPrio(blossom) -
2410                           (*_blossom_data)[blossom].offset){
2411                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2412                                   (*_blossom_data)[blossom].offset);
2413                }
2414              }
2415            }
2416          } else {
2417
2418            typename std::map<int, Arc>::iterator it =
2419              (*_node_data)[vi].heap_index.find(tree);
2420
2421            if (it != (*_node_data)[vi].heap_index.end()) {
2422              (*_node_data)[vi].heap.erase(it->second);
2423              (*_node_data)[vi].heap_index.erase(it);
2424              if ((*_node_data)[vi].heap.empty()) {
2425                _blossom_set->increase(v, std::numeric_limits<Value>::max());
2426              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
2427                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
2428              }
2429
2430              if ((*_blossom_data)[vb].status == MATCHED) {
2431                if (_blossom_set->classPrio(vb) ==
2432                    std::numeric_limits<Value>::max()) {
2433                  _delta2->erase(vb);
2434                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
2435                           (*_blossom_data)[vb].offset) {
2436                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
2437                                   (*_blossom_data)[vb].offset);
2438                }
2439              }
2440            }
2441          }
2442        }
2443      }
2444    }
2445
2446    void oddToMatched(int blossom) {
2447      (*_blossom_data)[blossom].offset -= _delta_sum;
2448
2449      if (_blossom_set->classPrio(blossom) !=
2450          std::numeric_limits<Value>::max()) {
2451        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2452                       (*_blossom_data)[blossom].offset);
2453      }
2454
2455      if (!_blossom_set->trivial(blossom)) {
2456        _delta4->erase(blossom);
2457      }
2458    }
2459
2460    void oddToEven(int blossom, int tree) {
2461      if (!_blossom_set->trivial(blossom)) {
2462        _delta4->erase(blossom);
2463        (*_blossom_data)[blossom].pot -=
2464          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
2465      }
2466
2467      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2468           n != INVALID; ++n) {
2469        int ni = (*_node_index)[n];
2470
2471        _blossom_set->increase(n, std::numeric_limits<Value>::max());
2472
2473        (*_node_data)[ni].heap.clear();
2474        (*_node_data)[ni].heap_index.clear();
2475        (*_node_data)[ni].pot +=
2476          2 * _delta_sum - (*_blossom_data)[blossom].offset;
2477
2478        for (InArcIt e(_graph, n); e != INVALID; ++e) {
2479          Node v = _graph.source(e);
2480          int vb = _blossom_set->find(v);
2481          int vi = (*_node_index)[v];
2482
2483          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2484            dualScale * _weight[e];
2485
2486          if ((*_blossom_data)[vb].status == EVEN) {
2487            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2488              _delta3->push(e, rw / 2);
2489            }
2490          } else {
2491
2492            typename std::map<int, Arc>::iterator it =
2493              (*_node_data)[vi].heap_index.find(tree);
2494
2495            if (it != (*_node_data)[vi].heap_index.end()) {
2496              if ((*_node_data)[vi].heap[it->second] > rw) {
2497                (*_node_data)[vi].heap.replace(it->second, e);
2498                (*_node_data)[vi].heap.decrease(e, rw);
2499                it->second = e;
2500              }
2501            } else {
2502              (*_node_data)[vi].heap.push(e, rw);
2503              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2504            }
2505
2506            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2507              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2508
2509              if ((*_blossom_data)[vb].status == MATCHED) {
2510                if (_delta2->state(vb) != _delta2->IN_HEAP) {
2511                  _delta2->push(vb, _blossom_set->classPrio(vb) -
2512                               (*_blossom_data)[vb].offset);
2513                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2514                           (*_blossom_data)[vb].offset) {
2515                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2516                                   (*_blossom_data)[vb].offset);
2517                }
2518              }
2519            }
2520          }
2521        }
2522      }
2523      (*_blossom_data)[blossom].offset = 0;
2524    }
2525
2526    void alternatePath(int even, int tree) {
2527      int odd;
2528
2529      evenToMatched(even, tree);
2530      (*_blossom_data)[even].status = MATCHED;
2531
2532      while ((*_blossom_data)[even].pred != INVALID) {
2533        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
2534        (*_blossom_data)[odd].status = MATCHED;
2535        oddToMatched(odd);
2536        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2537
2538        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
2539        (*_blossom_data)[even].status = MATCHED;
2540        evenToMatched(even, tree);
2541        (*_blossom_data)[even].next =
2542          _graph.oppositeArc((*_blossom_data)[odd].pred);
2543      }
2544
2545    }
2546
2547    void destroyTree(int tree) {
2548      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2549        if ((*_blossom_data)[b].status == EVEN) {
2550          (*_blossom_data)[b].status = MATCHED;
2551          evenToMatched(b, tree);
2552        } else if ((*_blossom_data)[b].status == ODD) {
2553          (*_blossom_data)[b].status = MATCHED;
2554          oddToMatched(b);
2555        }
2556      }
2557      _tree_set->eraseClass(tree);
2558    }
2559
2560    void augmentOnEdge(const Edge& edge) {
2561
2562      int left = _blossom_set->find(_graph.u(edge));
2563      int right = _blossom_set->find(_graph.v(edge));
2564
2565      int left_tree = _tree_set->find(left);
2566      alternatePath(left, left_tree);
2567      destroyTree(left_tree);
2568
2569      int right_tree = _tree_set->find(right);
2570      alternatePath(right, right_tree);
2571      destroyTree(right_tree);
2572
2573      (*_blossom_data)[left].next = _graph.direct(edge, true);
2574      (*_blossom_data)[right].next = _graph.direct(edge, false);
2575    }
2576
2577    void extendOnArc(const Arc& arc) {
2578      int base = _blossom_set->find(_graph.target(arc));
2579      int tree = _tree_set->find(base);
2580
2581      int odd = _blossom_set->find(_graph.source(arc));
2582      _tree_set->insert(odd, tree);
2583      (*_blossom_data)[odd].status = ODD;
2584      matchedToOdd(odd);
2585      (*_blossom_data)[odd].pred = arc;
2586
2587      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
2588      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2589      _tree_set->insert(even, tree);
2590      (*_blossom_data)[even].status = EVEN;
2591      matchedToEven(even, tree);
2592    }
2593
2594    void shrinkOnEdge(const Edge& edge, int tree) {
2595      int nca = -1;
2596      std::vector<int> left_path, right_path;
2597
2598      {
2599        std::set<int> left_set, right_set;
2600        int left = _blossom_set->find(_graph.u(edge));
2601        left_path.push_back(left);
2602        left_set.insert(left);
2603
2604        int right = _blossom_set->find(_graph.v(edge));
2605        right_path.push_back(right);
2606        right_set.insert(right);
2607
2608        while (true) {
2609
2610          if ((*_blossom_data)[left].pred == INVALID) break;
2611
2612          left =
2613            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2614          left_path.push_back(left);
2615          left =
2616            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2617          left_path.push_back(left);
2618
2619          left_set.insert(left);
2620
2621          if (right_set.find(left) != right_set.end()) {
2622            nca = left;
2623            break;
2624          }
2625
2626          if ((*_blossom_data)[right].pred == INVALID) break;
2627
2628          right =
2629            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2630          right_path.push_back(right);
2631          right =
2632            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2633          right_path.push_back(right);
2634
2635          right_set.insert(right);
2636
2637          if (left_set.find(right) != left_set.end()) {
2638            nca = right;
2639            break;
2640          }
2641
2642        }
2643
2644        if (nca == -1) {
2645          if ((*_blossom_data)[left].pred == INVALID) {
2646            nca = right;
2647            while (left_set.find(nca) == left_set.end()) {
2648              nca =
2649                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2650              right_path.push_back(nca);
2651              nca =
2652                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2653              right_path.push_back(nca);
2654            }
2655          } else {
2656            nca = left;
2657            while (right_set.find(nca) == right_set.end()) {
2658              nca =
2659                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2660              left_path.push_back(nca);
2661              nca =
2662                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2663              left_path.push_back(nca);
2664            }
2665          }
2666        }
2667      }
2668
2669      std::vector<int> subblossoms;
2670      Arc prev;
2671
2672      prev = _graph.direct(edge, true);
2673      for (int i = 0; left_path[i] != nca; i += 2) {
2674        subblossoms.push_back(left_path[i]);
2675        (*_blossom_data)[left_path[i]].next = prev;
2676        _tree_set->erase(left_path[i]);
2677
2678        subblossoms.push_back(left_path[i + 1]);
2679        (*_blossom_data)[left_path[i + 1]].status = EVEN;
2680        oddToEven(left_path[i + 1], tree);
2681        _tree_set->erase(left_path[i + 1]);
2682        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
2683      }
2684
2685      int k = 0;
2686      while (right_path[k] != nca) ++k;
2687
2688      subblossoms.push_back(nca);
2689      (*_blossom_data)[nca].next = prev;
2690
2691      for (int i = k - 2; i >= 0; i -= 2) {
2692        subblossoms.push_back(right_path[i + 1]);
2693        (*_blossom_data)[right_path[i + 1]].status = EVEN;
2694        oddToEven(right_path[i + 1], tree);
2695        _tree_set->erase(right_path[i + 1]);
2696
2697        (*_blossom_data)[right_path[i + 1]].next =
2698          (*_blossom_data)[right_path[i + 1]].pred;
2699
2700        subblossoms.push_back(right_path[i]);
2701        _tree_set->erase(right_path[i]);
2702      }
2703
2704      int surface =
2705        _blossom_set->join(subblossoms.begin(), subblossoms.end());
2706
2707      for (int i = 0; i < int(subblossoms.size()); ++i) {
2708        if (!_blossom_set->trivial(subblossoms[i])) {
2709          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
2710        }
2711        (*_blossom_data)[subblossoms[i]].status = MATCHED;
2712      }
2713
2714      (*_blossom_data)[surface].pot = -2 * _delta_sum;
2715      (*_blossom_data)[surface].offset = 0;
2716      (*_blossom_data)[surface].status = EVEN;
2717      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
2718      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
2719
2720      _tree_set->insert(surface, tree);
2721      _tree_set->erase(nca);
2722    }
2723
2724    void splitBlossom(int blossom) {
2725      Arc next = (*_blossom_data)[blossom].next;
2726      Arc pred = (*_blossom_data)[blossom].pred;
2727
2728      int tree = _tree_set->find(blossom);
2729
2730      (*_blossom_data)[blossom].status = MATCHED;
2731      oddToMatched(blossom);
2732      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2733        _delta2->erase(blossom);
2734      }
2735
2736      std::vector<int> subblossoms;
2737      _blossom_set->split(blossom, std::back_inserter(subblossoms));
2738
2739      Value offset = (*_blossom_data)[blossom].offset;
2740      int b = _blossom_set->find(_graph.source(pred));
2741      int d = _blossom_set->find(_graph.source(next));
2742
2743      int ib = -1, id = -1;
2744      for (int i = 0; i < int(subblossoms.size()); ++i) {
2745        if (subblossoms[i] == b) ib = i;
2746        if (subblossoms[i] == d) id = i;
2747
2748        (*_blossom_data)[subblossoms[i]].offset = offset;
2749        if (!_blossom_set->trivial(subblossoms[i])) {
2750          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2751        }
2752        if (_blossom_set->classPrio(subblossoms[i]) !=
2753            std::numeric_limits<Value>::max()) {
2754          _delta2->push(subblossoms[i],
2755                        _blossom_set->classPrio(subblossoms[i]) -
2756                        (*_blossom_data)[subblossoms[i]].offset);
2757        }
2758      }
2759
2760      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2761        for (int i = (id + 1) % subblossoms.size();
2762             i != ib; i = (i + 2) % subblossoms.size()) {
2763          int sb = subblossoms[i];
2764          int tb = subblossoms[(i + 1) % subblossoms.size()];
2765          (*_blossom_data)[sb].next =
2766            _graph.oppositeArc((*_blossom_data)[tb].next);
2767        }
2768
2769        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2770          int sb = subblossoms[i];
2771          int tb = subblossoms[(i + 1) % subblossoms.size()];
2772          int ub = subblossoms[(i + 2) % subblossoms.size()];
2773
2774          (*_blossom_data)[sb].status = ODD;
2775          matchedToOdd(sb);
2776          _tree_set->insert(sb, tree);
2777          (*_blossom_data)[sb].pred = pred;
2778          (*_blossom_data)[sb].next =
2779                           _graph.oppositeArc((*_blossom_data)[tb].next);
2780
2781          pred = (*_blossom_data)[ub].next;
2782
2783          (*_blossom_data)[tb].status = EVEN;
2784          matchedToEven(tb, tree);
2785          _tree_set->insert(tb, tree);
2786          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2787        }
2788
2789        (*_blossom_data)[subblossoms[id]].status = ODD;
2790        matchedToOdd(subblossoms[id]);
2791        _tree_set->insert(subblossoms[id], tree);
2792        (*_blossom_data)[subblossoms[id]].next = next;
2793        (*_blossom_data)[subblossoms[id]].pred = pred;
2794
2795      } else {
2796
2797        for (int i = (ib + 1) % subblossoms.size();
2798             i != id; i = (i + 2) % subblossoms.size()) {
2799          int sb = subblossoms[i];
2800          int tb = subblossoms[(i + 1) % subblossoms.size()];
2801          (*_blossom_data)[sb].next =
2802            _graph.oppositeArc((*_blossom_data)[tb].next);
2803        }
2804
2805        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2806          int sb = subblossoms[i];
2807          int tb = subblossoms[(i + 1) % subblossoms.size()];
2808          int ub = subblossoms[(i + 2) % subblossoms.size()];
2809
2810          (*_blossom_data)[sb].status = ODD;
2811          matchedToOdd(sb);
2812          _tree_set->insert(sb, tree);
2813          (*_blossom_data)[sb].next = next;
2814          (*_blossom_data)[sb].pred =
2815            _graph.oppositeArc((*_blossom_data)[tb].next);
2816
2817          (*_blossom_data)[tb].status = EVEN;
2818          matchedToEven(tb, tree);
2819          _tree_set->insert(tb, tree);
2820          (*_blossom_data)[tb].pred =
2821            (*_blossom_data)[tb].next =
2822            _graph.oppositeArc((*_blossom_data)[ub].next);
2823          next = (*_blossom_data)[ub].next;
2824        }
2825
2826        (*_blossom_data)[subblossoms[ib]].status = ODD;
2827        matchedToOdd(subblossoms[ib]);
2828        _tree_set->insert(subblossoms[ib], tree);
2829        (*_blossom_data)[subblossoms[ib]].next = next;
2830        (*_blossom_data)[subblossoms[ib]].pred = pred;
2831      }
2832      _tree_set->erase(blossom);
2833    }
2834
2835    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2836      if (_blossom_set->trivial(blossom)) {
2837        int bi = (*_node_index)[base];
2838        Value pot = (*_node_data)[bi].pot;
2839
2840        (*_matching)[base] = matching;
2841        _blossom_node_list.push_back(base);
2842        (*_node_potential)[base] = pot;
2843      } else {
2844
2845        Value pot = (*_blossom_data)[blossom].pot;
2846        int bn = _blossom_node_list.size();
2847
2848        std::vector<int> subblossoms;
2849        _blossom_set->split(blossom, std::back_inserter(subblossoms));
2850        int b = _blossom_set->find(base);
2851        int ib = -1;
2852        for (int i = 0; i < int(subblossoms.size()); ++i) {
2853          if (subblossoms[i] == b) { ib = i; break; }
2854        }
2855
2856        for (int i = 1; i < int(subblossoms.size()); i += 2) {
2857          int sb = subblossoms[(ib + i) % subblossoms.size()];
2858          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2859
2860          Arc m = (*_blossom_data)[tb].next;
2861          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2862          extractBlossom(tb, _graph.source(m), m);
2863        }
2864        extractBlossom(subblossoms[ib], base, matching);
2865
2866        int en = _blossom_node_list.size();
2867
2868        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2869      }
2870    }
2871
2872    void extractMatching() {
2873      std::vector<int> blossoms;
2874      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2875        blossoms.push_back(c);
2876      }
2877
2878      for (int i = 0; i < int(blossoms.size()); ++i) {
2879
2880        Value offset = (*_blossom_data)[blossoms[i]].offset;
2881        (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2882        for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2883             n != INVALID; ++n) {
2884          (*_node_data)[(*_node_index)[n]].pot -= offset;
2885        }
2886
2887        Arc matching = (*_blossom_data)[blossoms[i]].next;
2888        Node base = _graph.source(matching);
2889        extractBlossom(blossoms[i], base, matching);
2890      }
2891    }
2892
2893  public:
2894
2895    /// \brief Constructor
2896    ///
2897    /// Constructor.
2898    MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2899      : _graph(graph), _weight(weight), _matching(0),
2900        _node_potential(0), _blossom_potential(), _blossom_node_list(),
2901        _node_num(0), _blossom_num(0),
2902
2903        _blossom_index(0), _blossom_set(0), _blossom_data(0),
2904        _node_index(0), _node_heap_index(0), _node_data(0),
2905        _tree_set_index(0), _tree_set(0),
2906
2907        _delta2_index(0), _delta2(0),
2908        _delta3_index(0), _delta3(0),
2909        _delta4_index(0), _delta4(0),
2910
2911        _delta_sum() {}
2912
2913    ~MaxWeightedPerfectMatching() {
2914      destroyStructures();
2915    }
2916
2917    /// \name Execution Control
2918    /// The simplest way to execute the algorithm is to use the
2919    /// \ref run() member function.
2920
2921    ///@{
2922
2923    /// \brief Initialize the algorithm
2924    ///
2925    /// This function initializes the algorithm.
2926    void init() {
2927      createStructures();
2928
2929      for (ArcIt e(_graph); e != INVALID; ++e) {
2930        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
2931      }
2932      for (EdgeIt e(_graph); e != INVALID; ++e) {
2933        (*_delta3_index)[e] = _delta3->PRE_HEAP;
2934      }
2935      for (int i = 0; i < _blossom_num; ++i) {
2936        (*_delta2_index)[i] = _delta2->PRE_HEAP;
2937        (*_delta4_index)[i] = _delta4->PRE_HEAP;
2938      }
2939
2940      int index = 0;
2941      for (NodeIt n(_graph); n != INVALID; ++n) {
2942        Value max = - std::numeric_limits<Value>::max();
2943        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2944          if (_graph.target(e) == n) continue;
2945          if ((dualScale * _weight[e]) / 2 > max) {
2946            max = (dualScale * _weight[e]) / 2;
2947          }
2948        }
2949        (*_node_index)[n] = index;
2950        (*_node_data)[index].pot = max;
2951        int blossom =
2952          _blossom_set->insert(n, std::numeric_limits<Value>::max());
2953
2954        _tree_set->insert(blossom);
2955
2956        (*_blossom_data)[blossom].status = EVEN;
2957        (*_blossom_data)[blossom].pred = INVALID;
2958        (*_blossom_data)[blossom].next = INVALID;
2959        (*_blossom_data)[blossom].pot = 0;
2960        (*_blossom_data)[blossom].offset = 0;
2961        ++index;
2962      }
2963      for (EdgeIt e(_graph); e != INVALID; ++e) {
2964        int si = (*_node_index)[_graph.u(e)];
2965        int ti = (*_node_index)[_graph.v(e)];
2966        if (_graph.u(e) != _graph.v(e)) {
2967          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
2968                            dualScale * _weight[e]) / 2);
2969        }
2970      }
2971    }
2972
2973    /// \brief Start the algorithm
2974    ///
2975    /// This function starts the algorithm.
2976    ///
2977    /// \pre \ref init() must be called before using this function.
2978    bool start() {
2979      enum OpType {
2980        D2, D3, D4
2981      };
2982
2983      int unmatched = _node_num;
2984      while (unmatched > 0) {
2985        Value d2 = !_delta2->empty() ?
2986          _delta2->prio() : std::numeric_limits<Value>::max();
2987
2988        Value d3 = !_delta3->empty() ?
2989          _delta3->prio() : std::numeric_limits<Value>::max();
2990
2991        Value d4 = !_delta4->empty() ?
2992          _delta4->prio() : std::numeric_limits<Value>::max();
2993
2994        _delta_sum = d2; OpType ot = D2;
2995        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
2996        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
2997
2998        if (_delta_sum == std::numeric_limits<Value>::max()) {
2999          return false;
3000        }
3001
3002        switch (ot) {
3003        case D2:
3004          {
3005            int blossom = _delta2->top();
3006            Node n = _blossom_set->classTop(blossom);
3007            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
3008            extendOnArc(e);
3009          }
3010          break;
3011        case D3:
3012          {
3013            Edge e = _delta3->top();
3014
3015            int left_blossom = _blossom_set->find(_graph.u(e));
3016            int right_blossom = _blossom_set->find(_graph.v(e));
3017
3018            if (left_blossom == right_blossom) {
3019              _delta3->pop();
3020            } else {
3021              int left_tree = _tree_set->find(left_blossom);
3022              int right_tree = _tree_set->find(right_blossom);
3023
3024              if (left_tree == right_tree) {
3025                shrinkOnEdge(e, left_tree);
3026              } else {
3027                augmentOnEdge(e);
3028                unmatched -= 2;
3029              }
3030            }
3031          } break;
3032        case D4:
3033          splitBlossom(_delta4->top());
3034          break;
3035        }
3036      }
3037      extractMatching();
3038      return true;
3039    }
3040
3041    /// \brief Run the algorithm.
3042    ///
3043    /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
3044    ///
3045    /// \note mwpm.run() is just a shortcut of the following code.
3046    /// \code
3047    ///   mwpm.init();
3048    ///   mwpm.start();
3049    /// \endcode
3050    bool run() {
3051      init();
3052      return start();
3053    }
3054
3055    /// @}
3056
3057    /// \name Primal Solution
3058    /// Functions to get the primal solution, i.e. the maximum weighted
3059    /// perfect matching.\n
3060    /// Either \ref run() or \ref start() function should be called before
3061    /// using them.
3062
3063    /// @{
3064
3065    /// \brief Return the weight of the matching.
3066    ///
3067    /// This function returns the weight of the found matching.
3068    ///
3069    /// \pre Either run() or start() must be called before using this function.
3070    Value matchingWeight() const {
3071      Value sum = 0;
3072      for (NodeIt n(_graph); n != INVALID; ++n) {
3073        if ((*_matching)[n] != INVALID) {
3074          sum += _weight[(*_matching)[n]];
3075        }
3076      }
3077      return sum /= 2;
3078    }
3079
3080    /// \brief Return \c true if the given edge is in the matching.
3081    ///
3082    /// This function returns \c true if the given edge is in the found
3083    /// matching.
3084    ///
3085    /// \pre Either run() or start() must be called before using this function.
3086    bool matching(const Edge& edge) const {
3087      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
3088    }
3089
3090    /// \brief Return the matching arc (or edge) incident to the given node.
3091    ///
3092    /// This function returns the matching arc (or edge) incident to the
3093    /// given node in the found matching or \c INVALID if the node is
3094    /// not covered by the matching.
3095    ///
3096    /// \pre Either run() or start() must be called before using this function.
3097    Arc matching(const Node& node) const {
3098      return (*_matching)[node];
3099    }
3100
3101    /// \brief Return a const reference to the matching map.
3102    ///
3103    /// This function returns a const reference to a node map that stores
3104    /// the matching arc (or edge) incident to each node.
3105    const MatchingMap& matchingMap() const {
3106      return *_matching;
3107    }
3108
3109    /// \brief Return the mate of the given node.
3110    ///
3111    /// This function returns the mate of the given node in the found
3112    /// matching or \c INVALID if the node is not covered by the matching.
3113    ///
3114    /// \pre Either run() or start() must be called before using this function.
3115    Node mate(const Node& node) const {
3116      return _graph.target((*_matching)[node]);
3117    }
3118
3119    /// @}
3120
3121    /// \name Dual Solution
3122    /// Functions to get the dual solution.\n
3123    /// Either \ref run() or \ref start() function should be called before
3124    /// using them.
3125
3126    /// @{
3127
3128    /// \brief Return the value of the dual solution.
3129    ///
3130    /// This function returns the value of the dual solution.
3131    /// It should be equal to the primal value scaled by \ref dualScale
3132    /// "dual scale".
3133    ///
3134    /// \pre Either run() or start() must be called before using this function.
3135    Value dualValue() const {
3136      Value sum = 0;
3137      for (NodeIt n(_graph); n != INVALID; ++n) {
3138        sum += nodeValue(n);
3139      }
3140      for (int i = 0; i < blossomNum(); ++i) {
3141        sum += blossomValue(i) * (blossomSize(i) / 2);
3142      }
3143      return sum;
3144    }
3145
3146    /// \brief Return the dual value (potential) of the given node.
3147    ///
3148    /// This function returns the dual value (potential) of the given node.
3149    ///
3150    /// \pre Either run() or start() must be called before using this function.
3151    Value nodeValue(const Node& n) const {
3152      return (*_node_potential)[n];
3153    }
3154
3155    /// \brief Return the number of the blossoms in the basis.
3156    ///
3157    /// This function returns the number of the blossoms in the basis.
3158    ///
3159    /// \pre Either run() or start() must be called before using this function.
3160    /// \see BlossomIt
3161    int blossomNum() const {
3162      return _blossom_potential.size();
3163    }
3164
3165    /// \brief Return the number of the nodes in the given blossom.
3166    ///
3167    /// This function returns the number of the nodes in the given blossom.
3168    ///
3169    /// \pre Either run() or start() must be called before using this function.
3170    /// \see BlossomIt
3171    int blossomSize(int k) const {
3172      return _blossom_potential[k].end - _blossom_potential[k].begin;
3173    }
3174
3175    /// \brief Return the dual value (ptential) of the given blossom.
3176    ///
3177    /// This function returns the dual value (ptential) of the given blossom.
3178    ///
3179    /// \pre Either run() or start() must be called before using this function.
3180    Value blossomValue(int k) const {
3181      return _blossom_potential[k].value;
3182    }
3183
3184    /// \brief Iterator for obtaining the nodes of a blossom.
3185    ///
3186    /// This class provides an iterator for obtaining the nodes of the
3187    /// given blossom. It lists a subset of the nodes.
3188    /// Before using this iterator, you must allocate a
3189    /// MaxWeightedPerfectMatching class and execute it.
3190    class BlossomIt {
3191    public:
3192
3193      /// \brief Constructor.
3194      ///
3195      /// Constructor to get the nodes of the given variable.
3196      ///
3197      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
3198      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
3199      /// must be called before initializing this iterator.
3200      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
3201        : _algorithm(&algorithm)
3202      {
3203        _index = _algorithm->_blossom_potential[variable].begin;
3204        _last = _algorithm->_blossom_potential[variable].end;
3205      }
3206
3207      /// \brief Conversion to \c Node.
3208      ///
3209      /// Conversion to \c Node.
3210      operator Node() const {
3211        return _algorithm->_blossom_node_list[_index];
3212      }
3213
3214      /// \brief Increment operator.
3215      ///
3216      /// Increment operator.
3217      BlossomIt& operator++() {
3218        ++_index;
3219        return *this;
3220      }
3221
3222      /// \brief Validity checking
3223      ///
3224      /// This function checks whether the iterator is invalid.
3225      bool operator==(Invalid) const { return _index == _last; }
3226
3227      /// \brief Validity checking
3228      ///
3229      /// This function checks whether the iterator is valid.
3230      bool operator!=(Invalid) const { return _index != _last; }
3231
3232    private:
3233      const MaxWeightedPerfectMatching* _algorithm;
3234      int _last;
3235      int _index;
3236    };
3237
3238    /// @}
3239
3240  };
3241
3242} //END OF NAMESPACE LEMON
3243
3244#endif //LEMON_MAX_MATCHING_H
Note: See TracBrowser for help on using the repository browser.