COIN-OR::LEMON - Graph Library

source: lemon/lemon/matching.h @ 1081:f1398882a928

1.1 r1.1.4
Last change on this file since 1081:f1398882a928 was 1081:f1398882a928, checked in by Alpar Juttner <alpar@…>, 8 years ago

Unify sources

File size: 103.8 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-2011
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
809      if (!_node_potential) {
810        _node_potential = new NodePotential(_graph);
811      }
812
813      if (!_blossom_set) {
814        _blossom_index = new IntNodeMap(_graph);
815        _blossom_set = new BlossomSet(*_blossom_index);
816        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
817      } else if (_blossom_data->size() != _blossom_num) {
818        delete _blossom_data;
819        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
820      }
821
822      if (!_node_index) {
823        _node_index = new IntNodeMap(_graph);
824        _node_heap_index = new IntArcMap(_graph);
825        _node_data = new RangeMap<NodeData>(_node_num,
826                                            NodeData(*_node_heap_index));
827      } else {
828        delete _node_data;
829        _node_data = new RangeMap<NodeData>(_node_num,
830                                            NodeData(*_node_heap_index));
831      }
832
833      if (!_tree_set) {
834        _tree_set_index = new IntIntMap(_blossom_num);
835        _tree_set = new TreeSet(*_tree_set_index);
836      } else {
837        _tree_set_index->resize(_blossom_num);
838      }
839
840      if (!_delta1) {
841        _delta1_index = new IntNodeMap(_graph);
842        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
843      }
844
845      if (!_delta2) {
846        _delta2_index = new IntIntMap(_blossom_num);
847        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
848      } else {
849        _delta2_index->resize(_blossom_num);
850      }
851
852      if (!_delta3) {
853        _delta3_index = new IntEdgeMap(_graph);
854        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
855      }
856
857      if (!_delta4) {
858        _delta4_index = new IntIntMap(_blossom_num);
859        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
860      } else {
861        _delta4_index->resize(_blossom_num);
862      }
863    }
864
865    void destroyStructures() {
866      _node_num = countNodes(_graph);
867      _blossom_num = _node_num * 3 / 2;
868
869      if (_matching) {
870        delete _matching;
871      }
872      if (_node_potential) {
873        delete _node_potential;
874      }
875      if (_blossom_set) {
876        delete _blossom_index;
877        delete _blossom_set;
878        delete _blossom_data;
879      }
880
881      if (_node_index) {
882        delete _node_index;
883        delete _node_heap_index;
884        delete _node_data;
885      }
886
887      if (_tree_set) {
888        delete _tree_set_index;
889        delete _tree_set;
890      }
891      if (_delta1) {
892        delete _delta1_index;
893        delete _delta1;
894      }
895      if (_delta2) {
896        delete _delta2_index;
897        delete _delta2;
898      }
899      if (_delta3) {
900        delete _delta3_index;
901        delete _delta3;
902      }
903      if (_delta4) {
904        delete _delta4_index;
905        delete _delta4;
906      }
907    }
908
909    void matchedToEven(int blossom, int tree) {
910      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
911        _delta2->erase(blossom);
912      }
913
914      if (!_blossom_set->trivial(blossom)) {
915        (*_blossom_data)[blossom].pot -=
916          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
917      }
918
919      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
920           n != INVALID; ++n) {
921
922        _blossom_set->increase(n, std::numeric_limits<Value>::max());
923        int ni = (*_node_index)[n];
924
925        (*_node_data)[ni].heap.clear();
926        (*_node_data)[ni].heap_index.clear();
927
928        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
929
930        _delta1->push(n, (*_node_data)[ni].pot);
931
932        for (InArcIt e(_graph, n); e != INVALID; ++e) {
933          Node v = _graph.source(e);
934          int vb = _blossom_set->find(v);
935          int vi = (*_node_index)[v];
936
937          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
938            dualScale * _weight[e];
939
940          if ((*_blossom_data)[vb].status == EVEN) {
941            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
942              _delta3->push(e, rw / 2);
943            }
944          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
945            if (_delta3->state(e) != _delta3->IN_HEAP) {
946              _delta3->push(e, rw);
947            }
948          } else {
949            typename std::map<int, Arc>::iterator it =
950              (*_node_data)[vi].heap_index.find(tree);
951
952            if (it != (*_node_data)[vi].heap_index.end()) {
953              if ((*_node_data)[vi].heap[it->second] > rw) {
954                (*_node_data)[vi].heap.replace(it->second, e);
955                (*_node_data)[vi].heap.decrease(e, rw);
956                it->second = e;
957              }
958            } else {
959              (*_node_data)[vi].heap.push(e, rw);
960              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
961            }
962
963            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
964              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
965
966              if ((*_blossom_data)[vb].status == MATCHED) {
967                if (_delta2->state(vb) != _delta2->IN_HEAP) {
968                  _delta2->push(vb, _blossom_set->classPrio(vb) -
969                               (*_blossom_data)[vb].offset);
970                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
971                           (*_blossom_data)[vb].offset){
972                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
973                                   (*_blossom_data)[vb].offset);
974                }
975              }
976            }
977          }
978        }
979      }
980      (*_blossom_data)[blossom].offset = 0;
981    }
982
983    void matchedToOdd(int blossom) {
984      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
985        _delta2->erase(blossom);
986      }
987      (*_blossom_data)[blossom].offset += _delta_sum;
988      if (!_blossom_set->trivial(blossom)) {
989        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
990                     (*_blossom_data)[blossom].offset);
991      }
992    }
993
994    void evenToMatched(int blossom, int tree) {
995      if (!_blossom_set->trivial(blossom)) {
996        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
997      }
998
999      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1000           n != INVALID; ++n) {
1001        int ni = (*_node_index)[n];
1002        (*_node_data)[ni].pot -= _delta_sum;
1003
1004        _delta1->erase(n);
1005
1006        for (InArcIt e(_graph, n); e != INVALID; ++e) {
1007          Node v = _graph.source(e);
1008          int vb = _blossom_set->find(v);
1009          int vi = (*_node_index)[v];
1010
1011          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1012            dualScale * _weight[e];
1013
1014          if (vb == blossom) {
1015            if (_delta3->state(e) == _delta3->IN_HEAP) {
1016              _delta3->erase(e);
1017            }
1018          } else if ((*_blossom_data)[vb].status == EVEN) {
1019
1020            if (_delta3->state(e) == _delta3->IN_HEAP) {
1021              _delta3->erase(e);
1022            }
1023
1024            int vt = _tree_set->find(vb);
1025
1026            if (vt != tree) {
1027
1028              Arc r = _graph.oppositeArc(e);
1029
1030              typename std::map<int, Arc>::iterator it =
1031                (*_node_data)[ni].heap_index.find(vt);
1032
1033              if (it != (*_node_data)[ni].heap_index.end()) {
1034                if ((*_node_data)[ni].heap[it->second] > rw) {
1035                  (*_node_data)[ni].heap.replace(it->second, r);
1036                  (*_node_data)[ni].heap.decrease(r, rw);
1037                  it->second = r;
1038                }
1039              } else {
1040                (*_node_data)[ni].heap.push(r, rw);
1041                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1042              }
1043
1044              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1045                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1046
1047                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1048                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1049                               (*_blossom_data)[blossom].offset);
1050                } else if ((*_delta2)[blossom] >
1051                           _blossom_set->classPrio(blossom) -
1052                           (*_blossom_data)[blossom].offset){
1053                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1054                                   (*_blossom_data)[blossom].offset);
1055                }
1056              }
1057            }
1058
1059          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1060            if (_delta3->state(e) == _delta3->IN_HEAP) {
1061              _delta3->erase(e);
1062            }
1063          } else {
1064
1065            typename std::map<int, Arc>::iterator it =
1066              (*_node_data)[vi].heap_index.find(tree);
1067
1068            if (it != (*_node_data)[vi].heap_index.end()) {
1069              (*_node_data)[vi].heap.erase(it->second);
1070              (*_node_data)[vi].heap_index.erase(it);
1071              if ((*_node_data)[vi].heap.empty()) {
1072                _blossom_set->increase(v, std::numeric_limits<Value>::max());
1073              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1074                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1075              }
1076
1077              if ((*_blossom_data)[vb].status == MATCHED) {
1078                if (_blossom_set->classPrio(vb) ==
1079                    std::numeric_limits<Value>::max()) {
1080                  _delta2->erase(vb);
1081                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1082                           (*_blossom_data)[vb].offset) {
1083                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
1084                                   (*_blossom_data)[vb].offset);
1085                }
1086              }
1087            }
1088          }
1089        }
1090      }
1091    }
1092
1093    void oddToMatched(int blossom) {
1094      (*_blossom_data)[blossom].offset -= _delta_sum;
1095
1096      if (_blossom_set->classPrio(blossom) !=
1097          std::numeric_limits<Value>::max()) {
1098        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1099                       (*_blossom_data)[blossom].offset);
1100      }
1101
1102      if (!_blossom_set->trivial(blossom)) {
1103        _delta4->erase(blossom);
1104      }
1105    }
1106
1107    void oddToEven(int blossom, int tree) {
1108      if (!_blossom_set->trivial(blossom)) {
1109        _delta4->erase(blossom);
1110        (*_blossom_data)[blossom].pot -=
1111          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1112      }
1113
1114      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1115           n != INVALID; ++n) {
1116        int ni = (*_node_index)[n];
1117
1118        _blossom_set->increase(n, std::numeric_limits<Value>::max());
1119
1120        (*_node_data)[ni].heap.clear();
1121        (*_node_data)[ni].heap_index.clear();
1122        (*_node_data)[ni].pot +=
1123          2 * _delta_sum - (*_blossom_data)[blossom].offset;
1124
1125        _delta1->push(n, (*_node_data)[ni].pot);
1126
1127        for (InArcIt e(_graph, n); e != INVALID; ++e) {
1128          Node v = _graph.source(e);
1129          int vb = _blossom_set->find(v);
1130          int vi = (*_node_index)[v];
1131
1132          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1133            dualScale * _weight[e];
1134
1135          if ((*_blossom_data)[vb].status == EVEN) {
1136            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1137              _delta3->push(e, rw / 2);
1138            }
1139          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1140            if (_delta3->state(e) != _delta3->IN_HEAP) {
1141              _delta3->push(e, rw);
1142            }
1143          } else {
1144
1145            typename std::map<int, Arc>::iterator it =
1146              (*_node_data)[vi].heap_index.find(tree);
1147
1148            if (it != (*_node_data)[vi].heap_index.end()) {
1149              if ((*_node_data)[vi].heap[it->second] > rw) {
1150                (*_node_data)[vi].heap.replace(it->second, e);
1151                (*_node_data)[vi].heap.decrease(e, rw);
1152                it->second = e;
1153              }
1154            } else {
1155              (*_node_data)[vi].heap.push(e, rw);
1156              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1157            }
1158
1159            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1160              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1161
1162              if ((*_blossom_data)[vb].status == MATCHED) {
1163                if (_delta2->state(vb) != _delta2->IN_HEAP) {
1164                  _delta2->push(vb, _blossom_set->classPrio(vb) -
1165                               (*_blossom_data)[vb].offset);
1166                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1167                           (*_blossom_data)[vb].offset) {
1168                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1169                                   (*_blossom_data)[vb].offset);
1170                }
1171              }
1172            }
1173          }
1174        }
1175      }
1176      (*_blossom_data)[blossom].offset = 0;
1177    }
1178
1179
1180    void matchedToUnmatched(int blossom) {
1181      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1182        _delta2->erase(blossom);
1183      }
1184
1185      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1186           n != INVALID; ++n) {
1187        int ni = (*_node_index)[n];
1188
1189        _blossom_set->increase(n, std::numeric_limits<Value>::max());
1190
1191        (*_node_data)[ni].heap.clear();
1192        (*_node_data)[ni].heap_index.clear();
1193
1194        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1195          Node v = _graph.target(e);
1196          int vb = _blossom_set->find(v);
1197          int vi = (*_node_index)[v];
1198
1199          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1200            dualScale * _weight[e];
1201
1202          if ((*_blossom_data)[vb].status == EVEN) {
1203            if (_delta3->state(e) != _delta3->IN_HEAP) {
1204              _delta3->push(e, rw);
1205            }
1206          }
1207        }
1208      }
1209    }
1210
1211    void unmatchedToMatched(int blossom) {
1212      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1213           n != INVALID; ++n) {
1214        int ni = (*_node_index)[n];
1215
1216        for (InArcIt e(_graph, n); e != INVALID; ++e) {
1217          Node v = _graph.source(e);
1218          int vb = _blossom_set->find(v);
1219          int vi = (*_node_index)[v];
1220
1221          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1222            dualScale * _weight[e];
1223
1224          if (vb == blossom) {
1225            if (_delta3->state(e) == _delta3->IN_HEAP) {
1226              _delta3->erase(e);
1227            }
1228          } else if ((*_blossom_data)[vb].status == EVEN) {
1229
1230            if (_delta3->state(e) == _delta3->IN_HEAP) {
1231              _delta3->erase(e);
1232            }
1233
1234            int vt = _tree_set->find(vb);
1235
1236            Arc r = _graph.oppositeArc(e);
1237
1238            typename std::map<int, Arc>::iterator it =
1239              (*_node_data)[ni].heap_index.find(vt);
1240
1241            if (it != (*_node_data)[ni].heap_index.end()) {
1242              if ((*_node_data)[ni].heap[it->second] > rw) {
1243                (*_node_data)[ni].heap.replace(it->second, r);
1244                (*_node_data)[ni].heap.decrease(r, rw);
1245                it->second = r;
1246              }
1247            } else {
1248              (*_node_data)[ni].heap.push(r, rw);
1249              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1250            }
1251
1252            if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1253              _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1254
1255              if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1256                _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1257                             (*_blossom_data)[blossom].offset);
1258              } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
1259                         (*_blossom_data)[blossom].offset){
1260                _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1261                                 (*_blossom_data)[blossom].offset);
1262              }
1263            }
1264
1265          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1266            if (_delta3->state(e) == _delta3->IN_HEAP) {
1267              _delta3->erase(e);
1268            }
1269          }
1270        }
1271      }
1272    }
1273
1274    void alternatePath(int even, int tree) {
1275      int odd;
1276
1277      evenToMatched(even, tree);
1278      (*_blossom_data)[even].status = MATCHED;
1279
1280      while ((*_blossom_data)[even].pred != INVALID) {
1281        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
1282        (*_blossom_data)[odd].status = MATCHED;
1283        oddToMatched(odd);
1284        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1285
1286        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
1287        (*_blossom_data)[even].status = MATCHED;
1288        evenToMatched(even, tree);
1289        (*_blossom_data)[even].next =
1290          _graph.oppositeArc((*_blossom_data)[odd].pred);
1291      }
1292
1293    }
1294
1295    void destroyTree(int tree) {
1296      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1297        if ((*_blossom_data)[b].status == EVEN) {
1298          (*_blossom_data)[b].status = MATCHED;
1299          evenToMatched(b, tree);
1300        } else if ((*_blossom_data)[b].status == ODD) {
1301          (*_blossom_data)[b].status = MATCHED;
1302          oddToMatched(b);
1303        }
1304      }
1305      _tree_set->eraseClass(tree);
1306    }
1307
1308
1309    void unmatchNode(const Node& node) {
1310      int blossom = _blossom_set->find(node);
1311      int tree = _tree_set->find(blossom);
1312
1313      alternatePath(blossom, tree);
1314      destroyTree(tree);
1315
1316      (*_blossom_data)[blossom].status = UNMATCHED;
1317      (*_blossom_data)[blossom].base = node;
1318      matchedToUnmatched(blossom);
1319    }
1320
1321
1322    void augmentOnEdge(const Edge& edge) {
1323
1324      int left = _blossom_set->find(_graph.u(edge));
1325      int right = _blossom_set->find(_graph.v(edge));
1326
1327      if ((*_blossom_data)[left].status == EVEN) {
1328        int left_tree = _tree_set->find(left);
1329        alternatePath(left, left_tree);
1330        destroyTree(left_tree);
1331      } else {
1332        (*_blossom_data)[left].status = MATCHED;
1333        unmatchedToMatched(left);
1334      }
1335
1336      if ((*_blossom_data)[right].status == EVEN) {
1337        int right_tree = _tree_set->find(right);
1338        alternatePath(right, right_tree);
1339        destroyTree(right_tree);
1340      } else {
1341        (*_blossom_data)[right].status = MATCHED;
1342        unmatchedToMatched(right);
1343      }
1344
1345      (*_blossom_data)[left].next = _graph.direct(edge, true);
1346      (*_blossom_data)[right].next = _graph.direct(edge, false);
1347    }
1348
1349    void extendOnArc(const Arc& arc) {
1350      int base = _blossom_set->find(_graph.target(arc));
1351      int tree = _tree_set->find(base);
1352
1353      int odd = _blossom_set->find(_graph.source(arc));
1354      _tree_set->insert(odd, tree);
1355      (*_blossom_data)[odd].status = ODD;
1356      matchedToOdd(odd);
1357      (*_blossom_data)[odd].pred = arc;
1358
1359      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
1360      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1361      _tree_set->insert(even, tree);
1362      (*_blossom_data)[even].status = EVEN;
1363      matchedToEven(even, tree);
1364    }
1365
1366    void shrinkOnEdge(const Edge& edge, int tree) {
1367      int nca = -1;
1368      std::vector<int> left_path, right_path;
1369
1370      {
1371        std::set<int> left_set, right_set;
1372        int left = _blossom_set->find(_graph.u(edge));
1373        left_path.push_back(left);
1374        left_set.insert(left);
1375
1376        int right = _blossom_set->find(_graph.v(edge));
1377        right_path.push_back(right);
1378        right_set.insert(right);
1379
1380        while (true) {
1381
1382          if ((*_blossom_data)[left].pred == INVALID) break;
1383
1384          left =
1385            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1386          left_path.push_back(left);
1387          left =
1388            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1389          left_path.push_back(left);
1390
1391          left_set.insert(left);
1392
1393          if (right_set.find(left) != right_set.end()) {
1394            nca = left;
1395            break;
1396          }
1397
1398          if ((*_blossom_data)[right].pred == INVALID) break;
1399
1400          right =
1401            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1402          right_path.push_back(right);
1403          right =
1404            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1405          right_path.push_back(right);
1406
1407          right_set.insert(right);
1408
1409          if (left_set.find(right) != left_set.end()) {
1410            nca = right;
1411            break;
1412          }
1413
1414        }
1415
1416        if (nca == -1) {
1417          if ((*_blossom_data)[left].pred == INVALID) {
1418            nca = right;
1419            while (left_set.find(nca) == left_set.end()) {
1420              nca =
1421                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1422              right_path.push_back(nca);
1423              nca =
1424                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1425              right_path.push_back(nca);
1426            }
1427          } else {
1428            nca = left;
1429            while (right_set.find(nca) == right_set.end()) {
1430              nca =
1431                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1432              left_path.push_back(nca);
1433              nca =
1434                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1435              left_path.push_back(nca);
1436            }
1437          }
1438        }
1439      }
1440
1441      std::vector<int> subblossoms;
1442      Arc prev;
1443
1444      prev = _graph.direct(edge, true);
1445      for (int i = 0; left_path[i] != nca; i += 2) {
1446        subblossoms.push_back(left_path[i]);
1447        (*_blossom_data)[left_path[i]].next = prev;
1448        _tree_set->erase(left_path[i]);
1449
1450        subblossoms.push_back(left_path[i + 1]);
1451        (*_blossom_data)[left_path[i + 1]].status = EVEN;
1452        oddToEven(left_path[i + 1], tree);
1453        _tree_set->erase(left_path[i + 1]);
1454        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
1455      }
1456
1457      int k = 0;
1458      while (right_path[k] != nca) ++k;
1459
1460      subblossoms.push_back(nca);
1461      (*_blossom_data)[nca].next = prev;
1462
1463      for (int i = k - 2; i >= 0; i -= 2) {
1464        subblossoms.push_back(right_path[i + 1]);
1465        (*_blossom_data)[right_path[i + 1]].status = EVEN;
1466        oddToEven(right_path[i + 1], tree);
1467        _tree_set->erase(right_path[i + 1]);
1468
1469        (*_blossom_data)[right_path[i + 1]].next =
1470          (*_blossom_data)[right_path[i + 1]].pred;
1471
1472        subblossoms.push_back(right_path[i]);
1473        _tree_set->erase(right_path[i]);
1474      }
1475
1476      int surface =
1477        _blossom_set->join(subblossoms.begin(), subblossoms.end());
1478
1479      for (int i = 0; i < int(subblossoms.size()); ++i) {
1480        if (!_blossom_set->trivial(subblossoms[i])) {
1481          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1482        }
1483        (*_blossom_data)[subblossoms[i]].status = MATCHED;
1484      }
1485
1486      (*_blossom_data)[surface].pot = -2 * _delta_sum;
1487      (*_blossom_data)[surface].offset = 0;
1488      (*_blossom_data)[surface].status = EVEN;
1489      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1490      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1491
1492      _tree_set->insert(surface, tree);
1493      _tree_set->erase(nca);
1494    }
1495
1496    void splitBlossom(int blossom) {
1497      Arc next = (*_blossom_data)[blossom].next;
1498      Arc pred = (*_blossom_data)[blossom].pred;
1499
1500      int tree = _tree_set->find(blossom);
1501
1502      (*_blossom_data)[blossom].status = MATCHED;
1503      oddToMatched(blossom);
1504      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1505        _delta2->erase(blossom);
1506      }
1507
1508      std::vector<int> subblossoms;
1509      _blossom_set->split(blossom, std::back_inserter(subblossoms));
1510
1511      Value offset = (*_blossom_data)[blossom].offset;
1512      int b = _blossom_set->find(_graph.source(pred));
1513      int d = _blossom_set->find(_graph.source(next));
1514
1515      int ib = -1, id = -1;
1516      for (int i = 0; i < int(subblossoms.size()); ++i) {
1517        if (subblossoms[i] == b) ib = i;
1518        if (subblossoms[i] == d) id = i;
1519
1520        (*_blossom_data)[subblossoms[i]].offset = offset;
1521        if (!_blossom_set->trivial(subblossoms[i])) {
1522          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1523        }
1524        if (_blossom_set->classPrio(subblossoms[i]) !=
1525            std::numeric_limits<Value>::max()) {
1526          _delta2->push(subblossoms[i],
1527                        _blossom_set->classPrio(subblossoms[i]) -
1528                        (*_blossom_data)[subblossoms[i]].offset);
1529        }
1530      }
1531
1532      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1533        for (int i = (id + 1) % subblossoms.size();
1534             i != ib; i = (i + 2) % subblossoms.size()) {
1535          int sb = subblossoms[i];
1536          int tb = subblossoms[(i + 1) % subblossoms.size()];
1537          (*_blossom_data)[sb].next =
1538            _graph.oppositeArc((*_blossom_data)[tb].next);
1539        }
1540
1541        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1542          int sb = subblossoms[i];
1543          int tb = subblossoms[(i + 1) % subblossoms.size()];
1544          int ub = subblossoms[(i + 2) % subblossoms.size()];
1545
1546          (*_blossom_data)[sb].status = ODD;
1547          matchedToOdd(sb);
1548          _tree_set->insert(sb, tree);
1549          (*_blossom_data)[sb].pred = pred;
1550          (*_blossom_data)[sb].next =
1551                           _graph.oppositeArc((*_blossom_data)[tb].next);
1552
1553          pred = (*_blossom_data)[ub].next;
1554
1555          (*_blossom_data)[tb].status = EVEN;
1556          matchedToEven(tb, tree);
1557          _tree_set->insert(tb, tree);
1558          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1559        }
1560
1561        (*_blossom_data)[subblossoms[id]].status = ODD;
1562        matchedToOdd(subblossoms[id]);
1563        _tree_set->insert(subblossoms[id], tree);
1564        (*_blossom_data)[subblossoms[id]].next = next;
1565        (*_blossom_data)[subblossoms[id]].pred = pred;
1566
1567      } else {
1568
1569        for (int i = (ib + 1) % subblossoms.size();
1570             i != id; i = (i + 2) % subblossoms.size()) {
1571          int sb = subblossoms[i];
1572          int tb = subblossoms[(i + 1) % subblossoms.size()];
1573          (*_blossom_data)[sb].next =
1574            _graph.oppositeArc((*_blossom_data)[tb].next);
1575        }
1576
1577        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1578          int sb = subblossoms[i];
1579          int tb = subblossoms[(i + 1) % subblossoms.size()];
1580          int ub = subblossoms[(i + 2) % subblossoms.size()];
1581
1582          (*_blossom_data)[sb].status = ODD;
1583          matchedToOdd(sb);
1584          _tree_set->insert(sb, tree);
1585          (*_blossom_data)[sb].next = next;
1586          (*_blossom_data)[sb].pred =
1587            _graph.oppositeArc((*_blossom_data)[tb].next);
1588
1589          (*_blossom_data)[tb].status = EVEN;
1590          matchedToEven(tb, tree);
1591          _tree_set->insert(tb, tree);
1592          (*_blossom_data)[tb].pred =
1593            (*_blossom_data)[tb].next =
1594            _graph.oppositeArc((*_blossom_data)[ub].next);
1595          next = (*_blossom_data)[ub].next;
1596        }
1597
1598        (*_blossom_data)[subblossoms[ib]].status = ODD;
1599        matchedToOdd(subblossoms[ib]);
1600        _tree_set->insert(subblossoms[ib], tree);
1601        (*_blossom_data)[subblossoms[ib]].next = next;
1602        (*_blossom_data)[subblossoms[ib]].pred = pred;
1603      }
1604      _tree_set->erase(blossom);
1605    }
1606
1607    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1608      if (_blossom_set->trivial(blossom)) {
1609        int bi = (*_node_index)[base];
1610        Value pot = (*_node_data)[bi].pot;
1611
1612        (*_matching)[base] = matching;
1613        _blossom_node_list.push_back(base);
1614        (*_node_potential)[base] = pot;
1615      } else {
1616
1617        Value pot = (*_blossom_data)[blossom].pot;
1618        int bn = _blossom_node_list.size();
1619
1620        std::vector<int> subblossoms;
1621        _blossom_set->split(blossom, std::back_inserter(subblossoms));
1622        int b = _blossom_set->find(base);
1623        int ib = -1;
1624        for (int i = 0; i < int(subblossoms.size()); ++i) {
1625          if (subblossoms[i] == b) { ib = i; break; }
1626        }
1627
1628        for (int i = 1; i < int(subblossoms.size()); i += 2) {
1629          int sb = subblossoms[(ib + i) % subblossoms.size()];
1630          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1631
1632          Arc m = (*_blossom_data)[tb].next;
1633          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1634          extractBlossom(tb, _graph.source(m), m);
1635        }
1636        extractBlossom(subblossoms[ib], base, matching);
1637
1638        int en = _blossom_node_list.size();
1639
1640        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1641      }
1642    }
1643
1644    void extractMatching() {
1645      std::vector<int> blossoms;
1646      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1647        blossoms.push_back(c);
1648      }
1649
1650      for (int i = 0; i < int(blossoms.size()); ++i) {
1651        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
1652
1653          Value offset = (*_blossom_data)[blossoms[i]].offset;
1654          (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1655          for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1656               n != INVALID; ++n) {
1657            (*_node_data)[(*_node_index)[n]].pot -= offset;
1658          }
1659
1660          Arc matching = (*_blossom_data)[blossoms[i]].next;
1661          Node base = _graph.source(matching);
1662          extractBlossom(blossoms[i], base, matching);
1663        } else {
1664          Node base = (*_blossom_data)[blossoms[i]].base;
1665          extractBlossom(blossoms[i], base, INVALID);
1666        }
1667      }
1668    }
1669
1670  public:
1671
1672    /// \brief Constructor
1673    ///
1674    /// Constructor.
1675    MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1676      : _graph(graph), _weight(weight), _matching(0),
1677        _node_potential(0), _blossom_potential(), _blossom_node_list(),
1678        _node_num(0), _blossom_num(0),
1679
1680        _blossom_index(0), _blossom_set(0), _blossom_data(0),
1681        _node_index(0), _node_heap_index(0), _node_data(0),
1682        _tree_set_index(0), _tree_set(0),
1683
1684        _delta1_index(0), _delta1(0),
1685        _delta2_index(0), _delta2(0),
1686        _delta3_index(0), _delta3(0),
1687        _delta4_index(0), _delta4(0),
1688
1689        _delta_sum() {}
1690
1691    ~MaxWeightedMatching() {
1692      destroyStructures();
1693    }
1694
1695    /// \name Execution Control
1696    /// The simplest way to execute the algorithm is to use the
1697    /// \ref run() member function.
1698
1699    ///@{
1700
1701    /// \brief Initialize the algorithm
1702    ///
1703    /// This function initializes the algorithm.
1704    void init() {
1705      createStructures();
1706
1707      _blossom_node_list.clear();
1708      _blossom_potential.clear();
1709
1710      for (ArcIt e(_graph); e != INVALID; ++e) {
1711        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1712      }
1713      for (NodeIt n(_graph); n != INVALID; ++n) {
1714        (*_delta1_index)[n] = _delta1->PRE_HEAP;
1715      }
1716      for (EdgeIt e(_graph); e != INVALID; ++e) {
1717        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1718      }
1719      for (int i = 0; i < _blossom_num; ++i) {
1720        (*_delta2_index)[i] = _delta2->PRE_HEAP;
1721        (*_delta4_index)[i] = _delta4->PRE_HEAP;
1722      }
1723
1724      _delta1->clear();
1725      _delta2->clear();
1726      _delta3->clear();
1727      _delta4->clear();
1728      _blossom_set->clear();
1729      _tree_set->clear();
1730
1731      int index = 0;
1732      for (NodeIt n(_graph); n != INVALID; ++n) {
1733        Value max = 0;
1734        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1735          if (_graph.target(e) == n) continue;
1736          if ((dualScale * _weight[e]) / 2 > max) {
1737            max = (dualScale * _weight[e]) / 2;
1738          }
1739        }
1740        (*_node_index)[n] = index;
1741        (*_node_data)[index].heap_index.clear();
1742        (*_node_data)[index].heap.clear();
1743        (*_node_data)[index].pot = max;
1744        _delta1->push(n, max);
1745        int blossom =
1746          _blossom_set->insert(n, std::numeric_limits<Value>::max());
1747
1748        _tree_set->insert(blossom);
1749
1750        (*_blossom_data)[blossom].status = EVEN;
1751        (*_blossom_data)[blossom].pred = INVALID;
1752        (*_blossom_data)[blossom].next = INVALID;
1753        (*_blossom_data)[blossom].pot = 0;
1754        (*_blossom_data)[blossom].offset = 0;
1755        ++index;
1756      }
1757      for (EdgeIt e(_graph); e != INVALID; ++e) {
1758        int si = (*_node_index)[_graph.u(e)];
1759        int ti = (*_node_index)[_graph.v(e)];
1760        if (_graph.u(e) != _graph.v(e)) {
1761          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1762                            dualScale * _weight[e]) / 2);
1763        }
1764      }
1765    }
1766
1767    /// \brief Start the algorithm
1768    ///
1769    /// This function starts the algorithm.
1770    ///
1771    /// \pre \ref init() must be called before using this function.
1772    void start() {
1773      enum OpType {
1774        D1, D2, D3, D4
1775      };
1776
1777      int unmatched = _node_num;
1778      while (unmatched > 0) {
1779        Value d1 = !_delta1->empty() ?
1780          _delta1->prio() : std::numeric_limits<Value>::max();
1781
1782        Value d2 = !_delta2->empty() ?
1783          _delta2->prio() : std::numeric_limits<Value>::max();
1784
1785        Value d3 = !_delta3->empty() ?
1786          _delta3->prio() : std::numeric_limits<Value>::max();
1787
1788        Value d4 = !_delta4->empty() ?
1789          _delta4->prio() : std::numeric_limits<Value>::max();
1790
1791        _delta_sum = d1; OpType ot = D1;
1792        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1793        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1794        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1795
1796
1797        switch (ot) {
1798        case D1:
1799          {
1800            Node n = _delta1->top();
1801            unmatchNode(n);
1802            --unmatched;
1803          }
1804          break;
1805        case D2:
1806          {
1807            int blossom = _delta2->top();
1808            Node n = _blossom_set->classTop(blossom);
1809            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
1810            extendOnArc(e);
1811          }
1812          break;
1813        case D3:
1814          {
1815            Edge e = _delta3->top();
1816
1817            int left_blossom = _blossom_set->find(_graph.u(e));
1818            int right_blossom = _blossom_set->find(_graph.v(e));
1819
1820            if (left_blossom == right_blossom) {
1821              _delta3->pop();
1822            } else {
1823              int left_tree;
1824              if ((*_blossom_data)[left_blossom].status == EVEN) {
1825                left_tree = _tree_set->find(left_blossom);
1826              } else {
1827                left_tree = -1;
1828                ++unmatched;
1829              }
1830              int right_tree;
1831              if ((*_blossom_data)[right_blossom].status == EVEN) {
1832                right_tree = _tree_set->find(right_blossom);
1833              } else {
1834                right_tree = -1;
1835                ++unmatched;
1836              }
1837
1838              if (left_tree == right_tree) {
1839                shrinkOnEdge(e, left_tree);
1840              } else {
1841                augmentOnEdge(e);
1842                unmatched -= 2;
1843              }
1844            }
1845          } break;
1846        case D4:
1847          splitBlossom(_delta4->top());
1848          break;
1849        }
1850      }
1851      extractMatching();
1852    }
1853
1854    /// \brief Run the algorithm.
1855    ///
1856    /// This method runs the \c %MaxWeightedMatching algorithm.
1857    ///
1858    /// \note mwm.run() is just a shortcut of the following code.
1859    /// \code
1860    ///   mwm.init();
1861    ///   mwm.start();
1862    /// \endcode
1863    void run() {
1864      init();
1865      start();
1866    }
1867
1868    /// @}
1869
1870    /// \name Primal Solution
1871    /// Functions to get the primal solution, i.e. the maximum weighted
1872    /// matching.\n
1873    /// Either \ref run() or \ref start() function should be called before
1874    /// using them.
1875
1876    /// @{
1877
1878    /// \brief Return the weight of the matching.
1879    ///
1880    /// This function returns the weight of the found matching.
1881    ///
1882    /// \pre Either run() or start() must be called before using this function.
1883    Value matchingWeight() const {
1884      Value sum = 0;
1885      for (NodeIt n(_graph); n != INVALID; ++n) {
1886        if ((*_matching)[n] != INVALID) {
1887          sum += _weight[(*_matching)[n]];
1888        }
1889      }
1890      return sum /= 2;
1891    }
1892
1893    /// \brief Return the size (cardinality) of the matching.
1894    ///
1895    /// This function returns the size (cardinality) of the found matching.
1896    ///
1897    /// \pre Either run() or start() must be called before using this function.
1898    int matchingSize() const {
1899      int num = 0;
1900      for (NodeIt n(_graph); n != INVALID; ++n) {
1901        if ((*_matching)[n] != INVALID) {
1902          ++num;
1903        }
1904      }
1905      return num /= 2;
1906    }
1907
1908    /// \brief Return \c true if the given edge is in the matching.
1909    ///
1910    /// This function returns \c true if the given edge is in the found
1911    /// matching.
1912    ///
1913    /// \pre Either run() or start() must be called before using this function.
1914    bool matching(const Edge& edge) const {
1915      return edge == (*_matching)[_graph.u(edge)];
1916    }
1917
1918    /// \brief Return the matching arc (or edge) incident to the given node.
1919    ///
1920    /// This function returns the matching arc (or edge) incident to the
1921    /// given node in the found matching or \c INVALID if the node is
1922    /// not covered by the matching.
1923    ///
1924    /// \pre Either run() or start() must be called before using this function.
1925    Arc matching(const Node& node) const {
1926      return (*_matching)[node];
1927    }
1928
1929    /// \brief Return a const reference to the matching map.
1930    ///
1931    /// This function returns a const reference to a node map that stores
1932    /// the matching arc (or edge) incident to each node.
1933    const MatchingMap& matchingMap() const {
1934      return *_matching;
1935    }
1936
1937    /// \brief Return the mate of the given node.
1938    ///
1939    /// This function returns the mate of the given node in the found
1940    /// matching or \c INVALID if the node is not covered by the matching.
1941    ///
1942    /// \pre Either run() or start() must be called before using this function.
1943    Node mate(const Node& node) const {
1944      return (*_matching)[node] != INVALID ?
1945        _graph.target((*_matching)[node]) : INVALID;
1946    }
1947
1948    /// @}
1949
1950    /// \name Dual Solution
1951    /// Functions to get the dual solution.\n
1952    /// Either \ref run() or \ref start() function should be called before
1953    /// using them.
1954
1955    /// @{
1956
1957    /// \brief Return the value of the dual solution.
1958    ///
1959    /// This function returns the value of the dual solution.
1960    /// It should be equal to the primal value scaled by \ref dualScale
1961    /// "dual scale".
1962    ///
1963    /// \pre Either run() or start() must be called before using this function.
1964    Value dualValue() const {
1965      Value sum = 0;
1966      for (NodeIt n(_graph); n != INVALID; ++n) {
1967        sum += nodeValue(n);
1968      }
1969      for (int i = 0; i < blossomNum(); ++i) {
1970        sum += blossomValue(i) * (blossomSize(i) / 2);
1971      }
1972      return sum;
1973    }
1974
1975    /// \brief Return the dual value (potential) of the given node.
1976    ///
1977    /// This function returns the dual value (potential) of the given node.
1978    ///
1979    /// \pre Either run() or start() must be called before using this function.
1980    Value nodeValue(const Node& n) const {
1981      return (*_node_potential)[n];
1982    }
1983
1984    /// \brief Return the number of the blossoms in the basis.
1985    ///
1986    /// This function returns the number of the blossoms in the basis.
1987    ///
1988    /// \pre Either run() or start() must be called before using this function.
1989    /// \see BlossomIt
1990    int blossomNum() const {
1991      return _blossom_potential.size();
1992    }
1993
1994    /// \brief Return the number of the nodes in the given blossom.
1995    ///
1996    /// This function returns the number of the nodes in the given blossom.
1997    ///
1998    /// \pre Either run() or start() must be called before using this function.
1999    /// \see BlossomIt
2000    int blossomSize(int k) const {
2001      return _blossom_potential[k].end - _blossom_potential[k].begin;
2002    }
2003
2004    /// \brief Return the dual value (ptential) of the given blossom.
2005    ///
2006    /// This function returns the dual value (ptential) of the given blossom.
2007    ///
2008    /// \pre Either run() or start() must be called before using this function.
2009    Value blossomValue(int k) const {
2010      return _blossom_potential[k].value;
2011    }
2012
2013    /// \brief Iterator for obtaining the nodes of a blossom.
2014    ///
2015    /// This class provides an iterator for obtaining the nodes of the
2016    /// given blossom. It lists a subset of the nodes.
2017    /// Before using this iterator, you must allocate a
2018    /// MaxWeightedMatching class and execute it.
2019    class BlossomIt {
2020    public:
2021
2022      /// \brief Constructor.
2023      ///
2024      /// Constructor to get the nodes of the given variable.
2025      ///
2026      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
2027      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
2028      /// called before initializing this iterator.
2029      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
2030        : _algorithm(&algorithm)
2031      {
2032        _index = _algorithm->_blossom_potential[variable].begin;
2033        _last = _algorithm->_blossom_potential[variable].end;
2034      }
2035
2036      /// \brief Conversion to \c Node.
2037      ///
2038      /// Conversion to \c Node.
2039      operator Node() const {
2040        return _algorithm->_blossom_node_list[_index];
2041      }
2042
2043      /// \brief Increment operator.
2044      ///
2045      /// Increment operator.
2046      BlossomIt& operator++() {
2047        ++_index;
2048        return *this;
2049      }
2050
2051      /// \brief Validity checking
2052      ///
2053      /// Checks whether the iterator is invalid.
2054      bool operator==(Invalid) const { return _index == _last; }
2055
2056      /// \brief Validity checking
2057      ///
2058      /// Checks whether the iterator is valid.
2059      bool operator!=(Invalid) const { return _index != _last; }
2060
2061    private:
2062      const MaxWeightedMatching* _algorithm;
2063      int _last;
2064      int _index;
2065    };
2066
2067    /// @}
2068
2069  };
2070
2071  /// \ingroup matching
2072  ///
2073  /// \brief Weighted perfect matching in general graphs
2074  ///
2075  /// This class provides an efficient implementation of Edmond's
2076  /// maximum weighted perfect matching algorithm. The implementation
2077  /// is based on extensive use of priority queues and provides
2078  /// \f$O(nm\log n)\f$ time complexity.
2079  ///
2080  /// The maximum weighted perfect matching problem is to find a subset of
2081  /// the edges in an undirected graph with maximum overall weight for which
2082  /// each node has exactly one incident edge.
2083  /// It can be formulated with the following linear program.
2084  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
2085  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
2086      \quad \forall B\in\mathcal{O}\f] */
2087  /// \f[x_e \ge 0\quad \forall e\in E\f]
2088  /// \f[\max \sum_{e\in E}x_ew_e\f]
2089  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
2090  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
2091  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
2092  /// subsets of the nodes.
2093  ///
2094  /// The algorithm calculates an optimal matching and a proof of the
2095  /// optimality. The solution of the dual problem can be used to check
2096  /// the result of the algorithm. The dual linear problem is the
2097  /// following.
2098  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
2099      w_{uv} \quad \forall uv\in E\f] */
2100  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
2101  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
2102      \frac{\vert B \vert - 1}{2}z_B\f] */
2103  ///
2104  /// The algorithm can be executed with the run() function.
2105  /// After it the matching (the primal solution) and the dual solution
2106  /// can be obtained using the query functions and the
2107  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
2108  /// which is able to iterate on the nodes of a blossom.
2109  /// If the value type is integer, then the dual solution is multiplied
2110  /// by \ref MaxWeightedMatching::dualScale "4".
2111  ///
2112  /// \tparam GR The undirected graph type the algorithm runs on.
2113  /// \tparam WM The type edge weight map. The default type is
2114  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
2115#ifdef DOXYGEN
2116  template <typename GR, typename WM>
2117#else
2118  template <typename GR,
2119            typename WM = typename GR::template EdgeMap<int> >
2120#endif
2121  class MaxWeightedPerfectMatching {
2122  public:
2123
2124    /// The graph type of the algorithm
2125    typedef GR Graph;
2126    /// The type of the edge weight map
2127    typedef WM WeightMap;
2128    /// The value type of the edge weights
2129    typedef typename WeightMap::Value Value;
2130
2131    /// \brief Scaling factor for dual solution
2132    ///
2133    /// Scaling factor for dual solution, it is equal to 4 or 1
2134    /// according to the value type.
2135    static const int dualScale =
2136      std::numeric_limits<Value>::is_integer ? 4 : 1;
2137
2138    /// The type of the matching map
2139    typedef typename Graph::template NodeMap<typename Graph::Arc>
2140    MatchingMap;
2141
2142  private:
2143
2144    TEMPLATE_GRAPH_TYPEDEFS(Graph);
2145
2146    typedef typename Graph::template NodeMap<Value> NodePotential;
2147    typedef std::vector<Node> BlossomNodeList;
2148
2149    struct BlossomVariable {
2150      int begin, end;
2151      Value value;
2152
2153      BlossomVariable(int _begin, int _end, Value _value)
2154        : begin(_begin), end(_end), value(_value) {}
2155
2156    };
2157
2158    typedef std::vector<BlossomVariable> BlossomPotential;
2159
2160    const Graph& _graph;
2161    const WeightMap& _weight;
2162
2163    MatchingMap* _matching;
2164
2165    NodePotential* _node_potential;
2166
2167    BlossomPotential _blossom_potential;
2168    BlossomNodeList _blossom_node_list;
2169
2170    int _node_num;
2171    int _blossom_num;
2172
2173    typedef RangeMap<int> IntIntMap;
2174
2175    enum Status {
2176      EVEN = -1, MATCHED = 0, ODD = 1
2177    };
2178
2179    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
2180    struct BlossomData {
2181      int tree;
2182      Status status;
2183      Arc pred, next;
2184      Value pot, offset;
2185    };
2186
2187    IntNodeMap *_blossom_index;
2188    BlossomSet *_blossom_set;
2189    RangeMap<BlossomData>* _blossom_data;
2190
2191    IntNodeMap *_node_index;
2192    IntArcMap *_node_heap_index;
2193
2194    struct NodeData {
2195
2196      NodeData(IntArcMap& node_heap_index)
2197        : heap(node_heap_index) {}
2198
2199      int blossom;
2200      Value pot;
2201      BinHeap<Value, IntArcMap> heap;
2202      std::map<int, Arc> heap_index;
2203
2204      int tree;
2205    };
2206
2207    RangeMap<NodeData>* _node_data;
2208
2209    typedef ExtendFindEnum<IntIntMap> TreeSet;
2210
2211    IntIntMap *_tree_set_index;
2212    TreeSet *_tree_set;
2213
2214    IntIntMap *_delta2_index;
2215    BinHeap<Value, IntIntMap> *_delta2;
2216
2217    IntEdgeMap *_delta3_index;
2218    BinHeap<Value, IntEdgeMap> *_delta3;
2219
2220    IntIntMap *_delta4_index;
2221    BinHeap<Value, IntIntMap> *_delta4;
2222
2223    Value _delta_sum;
2224
2225    void createStructures() {
2226      _node_num = countNodes(_graph);
2227      _blossom_num = _node_num * 3 / 2;
2228
2229      if (!_matching) {
2230        _matching = new MatchingMap(_graph);
2231      }
2232
2233      if (!_node_potential) {
2234        _node_potential = new NodePotential(_graph);
2235      }
2236
2237      if (!_blossom_set) {
2238        _blossom_index = new IntNodeMap(_graph);
2239        _blossom_set = new BlossomSet(*_blossom_index);
2240        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2241      } else if (_blossom_data->size() != _blossom_num) {
2242        delete _blossom_data;
2243        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2244      }
2245
2246      if (!_node_index) {
2247        _node_index = new IntNodeMap(_graph);
2248        _node_heap_index = new IntArcMap(_graph);
2249        _node_data = new RangeMap<NodeData>(_node_num,
2250                                            NodeData(*_node_heap_index));
2251      } else if (_node_data->size() != _node_num) {
2252        delete _node_data;
2253        _node_data = new RangeMap<NodeData>(_node_num,
2254                                            NodeData(*_node_heap_index));
2255      }
2256
2257      if (!_tree_set) {
2258        _tree_set_index = new IntIntMap(_blossom_num);
2259        _tree_set = new TreeSet(*_tree_set_index);
2260      } else {
2261        _tree_set_index->resize(_blossom_num);
2262      }
2263
2264      if (!_delta2) {
2265        _delta2_index = new IntIntMap(_blossom_num);
2266        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2267      } else {
2268        _delta2_index->resize(_blossom_num);
2269      }
2270
2271      if (!_delta3) {
2272        _delta3_index = new IntEdgeMap(_graph);
2273        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
2274      }
2275
2276      if (!_delta4) {
2277        _delta4_index = new IntIntMap(_blossom_num);
2278        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2279      } else {
2280        _delta4_index->resize(_blossom_num);
2281      }
2282    }
2283
2284    void destroyStructures() {
2285      _node_num = countNodes(_graph);
2286      _blossom_num = _node_num * 3 / 2;
2287
2288      if (_matching) {
2289        delete _matching;
2290      }
2291      if (_node_potential) {
2292        delete _node_potential;
2293      }
2294      if (_blossom_set) {
2295        delete _blossom_index;
2296        delete _blossom_set;
2297        delete _blossom_data;
2298      }
2299
2300      if (_node_index) {
2301        delete _node_index;
2302        delete _node_heap_index;
2303        delete _node_data;
2304      }
2305
2306      if (_tree_set) {
2307        delete _tree_set_index;
2308        delete _tree_set;
2309      }
2310      if (_delta2) {
2311        delete _delta2_index;
2312        delete _delta2;
2313      }
2314      if (_delta3) {
2315        delete _delta3_index;
2316        delete _delta3;
2317      }
2318      if (_delta4) {
2319        delete _delta4_index;
2320        delete _delta4;
2321      }
2322    }
2323
2324    void matchedToEven(int blossom, int tree) {
2325      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2326        _delta2->erase(blossom);
2327      }
2328
2329      if (!_blossom_set->trivial(blossom)) {
2330        (*_blossom_data)[blossom].pot -=
2331          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2332      }
2333
2334      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2335           n != INVALID; ++n) {
2336
2337        _blossom_set->increase(n, std::numeric_limits<Value>::max());
2338        int ni = (*_node_index)[n];
2339
2340        (*_node_data)[ni].heap.clear();
2341        (*_node_data)[ni].heap_index.clear();
2342
2343        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
2344
2345        for (InArcIt e(_graph, n); e != INVALID; ++e) {
2346          Node v = _graph.source(e);
2347          int vb = _blossom_set->find(v);
2348          int vi = (*_node_index)[v];
2349
2350          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2351            dualScale * _weight[e];
2352
2353          if ((*_blossom_data)[vb].status == EVEN) {
2354            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2355              _delta3->push(e, rw / 2);
2356            }
2357          } else {
2358            typename std::map<int, Arc>::iterator it =
2359              (*_node_data)[vi].heap_index.find(tree);
2360
2361            if (it != (*_node_data)[vi].heap_index.end()) {
2362              if ((*_node_data)[vi].heap[it->second] > rw) {
2363                (*_node_data)[vi].heap.replace(it->second, e);
2364                (*_node_data)[vi].heap.decrease(e, rw);
2365                it->second = e;
2366              }
2367            } else {
2368              (*_node_data)[vi].heap.push(e, rw);
2369              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2370            }
2371
2372            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2373              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2374
2375              if ((*_blossom_data)[vb].status == MATCHED) {
2376                if (_delta2->state(vb) != _delta2->IN_HEAP) {
2377                  _delta2->push(vb, _blossom_set->classPrio(vb) -
2378                               (*_blossom_data)[vb].offset);
2379                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2380                           (*_blossom_data)[vb].offset){
2381                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2382                                   (*_blossom_data)[vb].offset);
2383                }
2384              }
2385            }
2386          }
2387        }
2388      }
2389      (*_blossom_data)[blossom].offset = 0;
2390    }
2391
2392    void matchedToOdd(int blossom) {
2393      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2394        _delta2->erase(blossom);
2395      }
2396      (*_blossom_data)[blossom].offset += _delta_sum;
2397      if (!_blossom_set->trivial(blossom)) {
2398        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
2399                     (*_blossom_data)[blossom].offset);
2400      }
2401    }
2402
2403    void evenToMatched(int blossom, int tree) {
2404      if (!_blossom_set->trivial(blossom)) {
2405        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
2406      }
2407
2408      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2409           n != INVALID; ++n) {
2410        int ni = (*_node_index)[n];
2411        (*_node_data)[ni].pot -= _delta_sum;
2412
2413        for (InArcIt e(_graph, n); e != INVALID; ++e) {
2414          Node v = _graph.source(e);
2415          int vb = _blossom_set->find(v);
2416          int vi = (*_node_index)[v];
2417
2418          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2419            dualScale * _weight[e];
2420
2421          if (vb == blossom) {
2422            if (_delta3->state(e) == _delta3->IN_HEAP) {
2423              _delta3->erase(e);
2424            }
2425          } else if ((*_blossom_data)[vb].status == EVEN) {
2426
2427            if (_delta3->state(e) == _delta3->IN_HEAP) {
2428              _delta3->erase(e);
2429            }
2430
2431            int vt = _tree_set->find(vb);
2432
2433            if (vt != tree) {
2434
2435              Arc r = _graph.oppositeArc(e);
2436
2437              typename std::map<int, Arc>::iterator it =
2438                (*_node_data)[ni].heap_index.find(vt);
2439
2440              if (it != (*_node_data)[ni].heap_index.end()) {
2441                if ((*_node_data)[ni].heap[it->second] > rw) {
2442                  (*_node_data)[ni].heap.replace(it->second, r);
2443                  (*_node_data)[ni].heap.decrease(r, rw);
2444                  it->second = r;
2445                }
2446              } else {
2447                (*_node_data)[ni].heap.push(r, rw);
2448                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2449              }
2450
2451              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2452                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2453
2454                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2455                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2456                               (*_blossom_data)[blossom].offset);
2457                } else if ((*_delta2)[blossom] >
2458                           _blossom_set->classPrio(blossom) -
2459                           (*_blossom_data)[blossom].offset){
2460                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2461                                   (*_blossom_data)[blossom].offset);
2462                }
2463              }
2464            }
2465          } else {
2466
2467            typename std::map<int, Arc>::iterator it =
2468              (*_node_data)[vi].heap_index.find(tree);
2469
2470            if (it != (*_node_data)[vi].heap_index.end()) {
2471              (*_node_data)[vi].heap.erase(it->second);
2472              (*_node_data)[vi].heap_index.erase(it);
2473              if ((*_node_data)[vi].heap.empty()) {
2474                _blossom_set->increase(v, std::numeric_limits<Value>::max());
2475              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
2476                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
2477              }
2478
2479              if ((*_blossom_data)[vb].status == MATCHED) {
2480                if (_blossom_set->classPrio(vb) ==
2481                    std::numeric_limits<Value>::max()) {
2482                  _delta2->erase(vb);
2483                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
2484                           (*_blossom_data)[vb].offset) {
2485                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
2486                                   (*_blossom_data)[vb].offset);
2487                }
2488              }
2489            }
2490          }
2491        }
2492      }
2493    }
2494
2495    void oddToMatched(int blossom) {
2496      (*_blossom_data)[blossom].offset -= _delta_sum;
2497
2498      if (_blossom_set->classPrio(blossom) !=
2499          std::numeric_limits<Value>::max()) {
2500        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2501                       (*_blossom_data)[blossom].offset);
2502      }
2503
2504      if (!_blossom_set->trivial(blossom)) {
2505        _delta4->erase(blossom);
2506      }
2507    }
2508
2509    void oddToEven(int blossom, int tree) {
2510      if (!_blossom_set->trivial(blossom)) {
2511        _delta4->erase(blossom);
2512        (*_blossom_data)[blossom].pot -=
2513          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
2514      }
2515
2516      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2517           n != INVALID; ++n) {
2518        int ni = (*_node_index)[n];
2519
2520        _blossom_set->increase(n, std::numeric_limits<Value>::max());
2521
2522        (*_node_data)[ni].heap.clear();
2523        (*_node_data)[ni].heap_index.clear();
2524        (*_node_data)[ni].pot +=
2525          2 * _delta_sum - (*_blossom_data)[blossom].offset;
2526
2527        for (InArcIt e(_graph, n); e != INVALID; ++e) {
2528          Node v = _graph.source(e);
2529          int vb = _blossom_set->find(v);
2530          int vi = (*_node_index)[v];
2531
2532          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2533            dualScale * _weight[e];
2534
2535          if ((*_blossom_data)[vb].status == EVEN) {
2536            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2537              _delta3->push(e, rw / 2);
2538            }
2539          } else {
2540
2541            typename std::map<int, Arc>::iterator it =
2542              (*_node_data)[vi].heap_index.find(tree);
2543
2544            if (it != (*_node_data)[vi].heap_index.end()) {
2545              if ((*_node_data)[vi].heap[it->second] > rw) {
2546                (*_node_data)[vi].heap.replace(it->second, e);
2547                (*_node_data)[vi].heap.decrease(e, rw);
2548                it->second = e;
2549              }
2550            } else {
2551              (*_node_data)[vi].heap.push(e, rw);
2552              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2553            }
2554
2555            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2556              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2557
2558              if ((*_blossom_data)[vb].status == MATCHED) {
2559                if (_delta2->state(vb) != _delta2->IN_HEAP) {
2560                  _delta2->push(vb, _blossom_set->classPrio(vb) -
2561                               (*_blossom_data)[vb].offset);
2562                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2563                           (*_blossom_data)[vb].offset) {
2564                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2565                                   (*_blossom_data)[vb].offset);
2566                }
2567              }
2568            }
2569          }
2570        }
2571      }
2572      (*_blossom_data)[blossom].offset = 0;
2573    }
2574
2575    void alternatePath(int even, int tree) {
2576      int odd;
2577
2578      evenToMatched(even, tree);
2579      (*_blossom_data)[even].status = MATCHED;
2580
2581      while ((*_blossom_data)[even].pred != INVALID) {
2582        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
2583        (*_blossom_data)[odd].status = MATCHED;
2584        oddToMatched(odd);
2585        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2586
2587        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
2588        (*_blossom_data)[even].status = MATCHED;
2589        evenToMatched(even, tree);
2590        (*_blossom_data)[even].next =
2591          _graph.oppositeArc((*_blossom_data)[odd].pred);
2592      }
2593
2594    }
2595
2596    void destroyTree(int tree) {
2597      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2598        if ((*_blossom_data)[b].status == EVEN) {
2599          (*_blossom_data)[b].status = MATCHED;
2600          evenToMatched(b, tree);
2601        } else if ((*_blossom_data)[b].status == ODD) {
2602          (*_blossom_data)[b].status = MATCHED;
2603          oddToMatched(b);
2604        }
2605      }
2606      _tree_set->eraseClass(tree);
2607    }
2608
2609    void augmentOnEdge(const Edge& edge) {
2610
2611      int left = _blossom_set->find(_graph.u(edge));
2612      int right = _blossom_set->find(_graph.v(edge));
2613
2614      int left_tree = _tree_set->find(left);
2615      alternatePath(left, left_tree);
2616      destroyTree(left_tree);
2617
2618      int right_tree = _tree_set->find(right);
2619      alternatePath(right, right_tree);
2620      destroyTree(right_tree);
2621
2622      (*_blossom_data)[left].next = _graph.direct(edge, true);
2623      (*_blossom_data)[right].next = _graph.direct(edge, false);
2624    }
2625
2626    void extendOnArc(const Arc& arc) {
2627      int base = _blossom_set->find(_graph.target(arc));
2628      int tree = _tree_set->find(base);
2629
2630      int odd = _blossom_set->find(_graph.source(arc));
2631      _tree_set->insert(odd, tree);
2632      (*_blossom_data)[odd].status = ODD;
2633      matchedToOdd(odd);
2634      (*_blossom_data)[odd].pred = arc;
2635
2636      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
2637      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2638      _tree_set->insert(even, tree);
2639      (*_blossom_data)[even].status = EVEN;
2640      matchedToEven(even, tree);
2641    }
2642
2643    void shrinkOnEdge(const Edge& edge, int tree) {
2644      int nca = -1;
2645      std::vector<int> left_path, right_path;
2646
2647      {
2648        std::set<int> left_set, right_set;
2649        int left = _blossom_set->find(_graph.u(edge));
2650        left_path.push_back(left);
2651        left_set.insert(left);
2652
2653        int right = _blossom_set->find(_graph.v(edge));
2654        right_path.push_back(right);
2655        right_set.insert(right);
2656
2657        while (true) {
2658
2659          if ((*_blossom_data)[left].pred == INVALID) break;
2660
2661          left =
2662            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2663          left_path.push_back(left);
2664          left =
2665            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2666          left_path.push_back(left);
2667
2668          left_set.insert(left);
2669
2670          if (right_set.find(left) != right_set.end()) {
2671            nca = left;
2672            break;
2673          }
2674
2675          if ((*_blossom_data)[right].pred == INVALID) break;
2676
2677          right =
2678            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2679          right_path.push_back(right);
2680          right =
2681            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2682          right_path.push_back(right);
2683
2684          right_set.insert(right);
2685
2686          if (left_set.find(right) != left_set.end()) {
2687            nca = right;
2688            break;
2689          }
2690
2691        }
2692
2693        if (nca == -1) {
2694          if ((*_blossom_data)[left].pred == INVALID) {
2695            nca = right;
2696            while (left_set.find(nca) == left_set.end()) {
2697              nca =
2698                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2699              right_path.push_back(nca);
2700              nca =
2701                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2702              right_path.push_back(nca);
2703            }
2704          } else {
2705            nca = left;
2706            while (right_set.find(nca) == right_set.end()) {
2707              nca =
2708                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2709              left_path.push_back(nca);
2710              nca =
2711                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2712              left_path.push_back(nca);
2713            }
2714          }
2715        }
2716      }
2717
2718      std::vector<int> subblossoms;
2719      Arc prev;
2720
2721      prev = _graph.direct(edge, true);
2722      for (int i = 0; left_path[i] != nca; i += 2) {
2723        subblossoms.push_back(left_path[i]);
2724        (*_blossom_data)[left_path[i]].next = prev;
2725        _tree_set->erase(left_path[i]);
2726
2727        subblossoms.push_back(left_path[i + 1]);
2728        (*_blossom_data)[left_path[i + 1]].status = EVEN;
2729        oddToEven(left_path[i + 1], tree);
2730        _tree_set->erase(left_path[i + 1]);
2731        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
2732      }
2733
2734      int k = 0;
2735      while (right_path[k] != nca) ++k;
2736
2737      subblossoms.push_back(nca);
2738      (*_blossom_data)[nca].next = prev;
2739
2740      for (int i = k - 2; i >= 0; i -= 2) {
2741        subblossoms.push_back(right_path[i + 1]);
2742        (*_blossom_data)[right_path[i + 1]].status = EVEN;
2743        oddToEven(right_path[i + 1], tree);
2744        _tree_set->erase(right_path[i + 1]);
2745
2746        (*_blossom_data)[right_path[i + 1]].next =
2747          (*_blossom_data)[right_path[i + 1]].pred;
2748
2749        subblossoms.push_back(right_path[i]);
2750        _tree_set->erase(right_path[i]);
2751      }
2752
2753      int surface =
2754        _blossom_set->join(subblossoms.begin(), subblossoms.end());
2755
2756      for (int i = 0; i < int(subblossoms.size()); ++i) {
2757        if (!_blossom_set->trivial(subblossoms[i])) {
2758          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
2759        }
2760        (*_blossom_data)[subblossoms[i]].status = MATCHED;
2761      }
2762
2763      (*_blossom_data)[surface].pot = -2 * _delta_sum;
2764      (*_blossom_data)[surface].offset = 0;
2765      (*_blossom_data)[surface].status = EVEN;
2766      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
2767      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
2768
2769      _tree_set->insert(surface, tree);
2770      _tree_set->erase(nca);
2771    }
2772
2773    void splitBlossom(int blossom) {
2774      Arc next = (*_blossom_data)[blossom].next;
2775      Arc pred = (*_blossom_data)[blossom].pred;
2776
2777      int tree = _tree_set->find(blossom);
2778
2779      (*_blossom_data)[blossom].status = MATCHED;
2780      oddToMatched(blossom);
2781      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2782        _delta2->erase(blossom);
2783      }
2784
2785      std::vector<int> subblossoms;
2786      _blossom_set->split(blossom, std::back_inserter(subblossoms));
2787
2788      Value offset = (*_blossom_data)[blossom].offset;
2789      int b = _blossom_set->find(_graph.source(pred));
2790      int d = _blossom_set->find(_graph.source(next));
2791
2792      int ib = -1, id = -1;
2793      for (int i = 0; i < int(subblossoms.size()); ++i) {
2794        if (subblossoms[i] == b) ib = i;
2795        if (subblossoms[i] == d) id = i;
2796
2797        (*_blossom_data)[subblossoms[i]].offset = offset;
2798        if (!_blossom_set->trivial(subblossoms[i])) {
2799          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2800        }
2801        if (_blossom_set->classPrio(subblossoms[i]) !=
2802            std::numeric_limits<Value>::max()) {
2803          _delta2->push(subblossoms[i],
2804                        _blossom_set->classPrio(subblossoms[i]) -
2805                        (*_blossom_data)[subblossoms[i]].offset);
2806        }
2807      }
2808
2809      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2810        for (int i = (id + 1) % subblossoms.size();
2811             i != ib; i = (i + 2) % subblossoms.size()) {
2812          int sb = subblossoms[i];
2813          int tb = subblossoms[(i + 1) % subblossoms.size()];
2814          (*_blossom_data)[sb].next =
2815            _graph.oppositeArc((*_blossom_data)[tb].next);
2816        }
2817
2818        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2819          int sb = subblossoms[i];
2820          int tb = subblossoms[(i + 1) % subblossoms.size()];
2821          int ub = subblossoms[(i + 2) % subblossoms.size()];
2822
2823          (*_blossom_data)[sb].status = ODD;
2824          matchedToOdd(sb);
2825          _tree_set->insert(sb, tree);
2826          (*_blossom_data)[sb].pred = pred;
2827          (*_blossom_data)[sb].next =
2828                           _graph.oppositeArc((*_blossom_data)[tb].next);
2829
2830          pred = (*_blossom_data)[ub].next;
2831
2832          (*_blossom_data)[tb].status = EVEN;
2833          matchedToEven(tb, tree);
2834          _tree_set->insert(tb, tree);
2835          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2836        }
2837
2838        (*_blossom_data)[subblossoms[id]].status = ODD;
2839        matchedToOdd(subblossoms[id]);
2840        _tree_set->insert(subblossoms[id], tree);
2841        (*_blossom_data)[subblossoms[id]].next = next;
2842        (*_blossom_data)[subblossoms[id]].pred = pred;
2843
2844      } else {
2845
2846        for (int i = (ib + 1) % subblossoms.size();
2847             i != id; i = (i + 2) % subblossoms.size()) {
2848          int sb = subblossoms[i];
2849          int tb = subblossoms[(i + 1) % subblossoms.size()];
2850          (*_blossom_data)[sb].next =
2851            _graph.oppositeArc((*_blossom_data)[tb].next);
2852        }
2853
2854        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2855          int sb = subblossoms[i];
2856          int tb = subblossoms[(i + 1) % subblossoms.size()];
2857          int ub = subblossoms[(i + 2) % subblossoms.size()];
2858
2859          (*_blossom_data)[sb].status = ODD;
2860          matchedToOdd(sb);
2861          _tree_set->insert(sb, tree);
2862          (*_blossom_data)[sb].next = next;
2863          (*_blossom_data)[sb].pred =
2864            _graph.oppositeArc((*_blossom_data)[tb].next);
2865
2866          (*_blossom_data)[tb].status = EVEN;
2867          matchedToEven(tb, tree);
2868          _tree_set->insert(tb, tree);
2869          (*_blossom_data)[tb].pred =
2870            (*_blossom_data)[tb].next =
2871            _graph.oppositeArc((*_blossom_data)[ub].next);
2872          next = (*_blossom_data)[ub].next;
2873        }
2874
2875        (*_blossom_data)[subblossoms[ib]].status = ODD;
2876        matchedToOdd(subblossoms[ib]);
2877        _tree_set->insert(subblossoms[ib], tree);
2878        (*_blossom_data)[subblossoms[ib]].next = next;
2879        (*_blossom_data)[subblossoms[ib]].pred = pred;
2880      }
2881      _tree_set->erase(blossom);
2882    }
2883
2884    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2885      if (_blossom_set->trivial(blossom)) {
2886        int bi = (*_node_index)[base];
2887        Value pot = (*_node_data)[bi].pot;
2888
2889        (*_matching)[base] = matching;
2890        _blossom_node_list.push_back(base);
2891        (*_node_potential)[base] = pot;
2892      } else {
2893
2894        Value pot = (*_blossom_data)[blossom].pot;
2895        int bn = _blossom_node_list.size();
2896
2897        std::vector<int> subblossoms;
2898        _blossom_set->split(blossom, std::back_inserter(subblossoms));
2899        int b = _blossom_set->find(base);
2900        int ib = -1;
2901        for (int i = 0; i < int(subblossoms.size()); ++i) {
2902          if (subblossoms[i] == b) { ib = i; break; }
2903        }
2904
2905        for (int i = 1; i < int(subblossoms.size()); i += 2) {
2906          int sb = subblossoms[(ib + i) % subblossoms.size()];
2907          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2908
2909          Arc m = (*_blossom_data)[tb].next;
2910          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2911          extractBlossom(tb, _graph.source(m), m);
2912        }
2913        extractBlossom(subblossoms[ib], base, matching);
2914
2915        int en = _blossom_node_list.size();
2916
2917        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2918      }
2919    }
2920
2921    void extractMatching() {
2922      std::vector<int> blossoms;
2923      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2924        blossoms.push_back(c);
2925      }
2926
2927      for (int i = 0; i < int(blossoms.size()); ++i) {
2928
2929        Value offset = (*_blossom_data)[blossoms[i]].offset;
2930        (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2931        for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2932             n != INVALID; ++n) {
2933          (*_node_data)[(*_node_index)[n]].pot -= offset;
2934        }
2935
2936        Arc matching = (*_blossom_data)[blossoms[i]].next;
2937        Node base = _graph.source(matching);
2938        extractBlossom(blossoms[i], base, matching);
2939      }
2940    }
2941
2942  public:
2943
2944    /// \brief Constructor
2945    ///
2946    /// Constructor.
2947    MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2948      : _graph(graph), _weight(weight), _matching(0),
2949        _node_potential(0), _blossom_potential(), _blossom_node_list(),
2950        _node_num(0), _blossom_num(0),
2951
2952        _blossom_index(0), _blossom_set(0), _blossom_data(0),
2953        _node_index(0), _node_heap_index(0), _node_data(0),
2954        _tree_set_index(0), _tree_set(0),
2955
2956        _delta2_index(0), _delta2(0),
2957        _delta3_index(0), _delta3(0),
2958        _delta4_index(0), _delta4(0),
2959
2960        _delta_sum() {}
2961
2962    ~MaxWeightedPerfectMatching() {
2963      destroyStructures();
2964    }
2965
2966    /// \name Execution Control
2967    /// The simplest way to execute the algorithm is to use the
2968    /// \ref run() member function.
2969
2970    ///@{
2971
2972    /// \brief Initialize the algorithm
2973    ///
2974    /// This function initializes the algorithm.
2975    void init() {
2976      createStructures();
2977
2978      _blossom_node_list.clear();
2979      _blossom_potential.clear();
2980
2981      for (ArcIt e(_graph); e != INVALID; ++e) {
2982        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
2983      }
2984      for (EdgeIt e(_graph); e != INVALID; ++e) {
2985        (*_delta3_index)[e] = _delta3->PRE_HEAP;
2986      }
2987      for (int i = 0; i < _blossom_num; ++i) {
2988        (*_delta2_index)[i] = _delta2->PRE_HEAP;
2989        (*_delta4_index)[i] = _delta4->PRE_HEAP;
2990      }
2991
2992      _delta2->clear();
2993      _delta3->clear();
2994      _delta4->clear();
2995      _blossom_set->clear();
2996      _tree_set->clear();
2997
2998      int index = 0;
2999      for (NodeIt n(_graph); n != INVALID; ++n) {
3000        Value max = - std::numeric_limits<Value>::max();
3001        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
3002          if (_graph.target(e) == n) continue;
3003          if ((dualScale * _weight[e]) / 2 > max) {
3004            max = (dualScale * _weight[e]) / 2;
3005          }
3006        }
3007        (*_node_index)[n] = index;
3008        (*_node_data)[index].heap_index.clear();
3009        (*_node_data)[index].heap.clear();
3010        (*_node_data)[index].pot = max;
3011        int blossom =
3012          _blossom_set->insert(n, std::numeric_limits<Value>::max());
3013
3014        _tree_set->insert(blossom);
3015
3016        (*_blossom_data)[blossom].status = EVEN;
3017        (*_blossom_data)[blossom].pred = INVALID;
3018        (*_blossom_data)[blossom].next = INVALID;
3019        (*_blossom_data)[blossom].pot = 0;
3020        (*_blossom_data)[blossom].offset = 0;
3021        ++index;
3022      }
3023      for (EdgeIt e(_graph); e != INVALID; ++e) {
3024        int si = (*_node_index)[_graph.u(e)];
3025        int ti = (*_node_index)[_graph.v(e)];
3026        if (_graph.u(e) != _graph.v(e)) {
3027          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
3028                            dualScale * _weight[e]) / 2);
3029        }
3030      }
3031    }
3032
3033    /// \brief Start the algorithm
3034    ///
3035    /// This function starts the algorithm.
3036    ///
3037    /// \pre \ref init() must be called before using this function.
3038    bool start() {
3039      enum OpType {
3040        D2, D3, D4
3041      };
3042
3043      int unmatched = _node_num;
3044      while (unmatched > 0) {
3045        Value d2 = !_delta2->empty() ?
3046          _delta2->prio() : std::numeric_limits<Value>::max();
3047
3048        Value d3 = !_delta3->empty() ?
3049          _delta3->prio() : std::numeric_limits<Value>::max();
3050
3051        Value d4 = !_delta4->empty() ?
3052          _delta4->prio() : std::numeric_limits<Value>::max();
3053
3054        _delta_sum = d2; OpType ot = D2;
3055        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
3056        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
3057
3058        if (_delta_sum == std::numeric_limits<Value>::max()) {
3059          return false;
3060        }
3061
3062        switch (ot) {
3063        case D2:
3064          {
3065            int blossom = _delta2->top();
3066            Node n = _blossom_set->classTop(blossom);
3067            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
3068            extendOnArc(e);
3069          }
3070          break;
3071        case D3:
3072          {
3073            Edge e = _delta3->top();
3074
3075            int left_blossom = _blossom_set->find(_graph.u(e));
3076            int right_blossom = _blossom_set->find(_graph.v(e));
3077
3078            if (left_blossom == right_blossom) {
3079              _delta3->pop();
3080            } else {
3081              int left_tree = _tree_set->find(left_blossom);
3082              int right_tree = _tree_set->find(right_blossom);
3083
3084              if (left_tree == right_tree) {
3085                shrinkOnEdge(e, left_tree);
3086              } else {
3087                augmentOnEdge(e);
3088                unmatched -= 2;
3089              }
3090            }
3091          } break;
3092        case D4:
3093          splitBlossom(_delta4->top());
3094          break;
3095        }
3096      }
3097      extractMatching();
3098      return true;
3099    }
3100
3101    /// \brief Run the algorithm.
3102    ///
3103    /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
3104    ///
3105    /// \note mwpm.run() is just a shortcut of the following code.
3106    /// \code
3107    ///   mwpm.init();
3108    ///   mwpm.start();
3109    /// \endcode
3110    bool run() {
3111      init();
3112      return start();
3113    }
3114
3115    /// @}
3116
3117    /// \name Primal Solution
3118    /// Functions to get the primal solution, i.e. the maximum weighted
3119    /// perfect matching.\n
3120    /// Either \ref run() or \ref start() function should be called before
3121    /// using them.
3122
3123    /// @{
3124
3125    /// \brief Return the weight of the matching.
3126    ///
3127    /// This function returns the weight of the found matching.
3128    ///
3129    /// \pre Either run() or start() must be called before using this function.
3130    Value matchingWeight() const {
3131      Value sum = 0;
3132      for (NodeIt n(_graph); n != INVALID; ++n) {
3133        if ((*_matching)[n] != INVALID) {
3134          sum += _weight[(*_matching)[n]];
3135        }
3136      }
3137      return sum /= 2;
3138    }
3139
3140    /// \brief Return \c true if the given edge is in the matching.
3141    ///
3142    /// This function returns \c true if the given edge is in the found
3143    /// matching.
3144    ///
3145    /// \pre Either run() or start() must be called before using this function.
3146    bool matching(const Edge& edge) const {
3147      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
3148    }
3149
3150    /// \brief Return the matching arc (or edge) incident to the given node.
3151    ///
3152    /// This function returns the matching arc (or edge) incident to the
3153    /// given node in the found matching or \c INVALID if the node is
3154    /// not covered by the matching.
3155    ///
3156    /// \pre Either run() or start() must be called before using this function.
3157    Arc matching(const Node& node) const {
3158      return (*_matching)[node];
3159    }
3160
3161    /// \brief Return a const reference to the matching map.
3162    ///
3163    /// This function returns a const reference to a node map that stores
3164    /// the matching arc (or edge) incident to each node.
3165    const MatchingMap& matchingMap() const {
3166      return *_matching;
3167    }
3168
3169    /// \brief Return the mate of the given node.
3170    ///
3171    /// This function returns the mate of the given node in the found
3172    /// matching or \c INVALID if the node is not covered by the matching.
3173    ///
3174    /// \pre Either run() or start() must be called before using this function.
3175    Node mate(const Node& node) const {
3176      return _graph.target((*_matching)[node]);
3177    }
3178
3179    /// @}
3180
3181    /// \name Dual Solution
3182    /// Functions to get the dual solution.\n
3183    /// Either \ref run() or \ref start() function should be called before
3184    /// using them.
3185
3186    /// @{
3187
3188    /// \brief Return the value of the dual solution.
3189    ///
3190    /// This function returns the value of the dual solution.
3191    /// It should be equal to the primal value scaled by \ref dualScale
3192    /// "dual scale".
3193    ///
3194    /// \pre Either run() or start() must be called before using this function.
3195    Value dualValue() const {
3196      Value sum = 0;
3197      for (NodeIt n(_graph); n != INVALID; ++n) {
3198        sum += nodeValue(n);
3199      }
3200      for (int i = 0; i < blossomNum(); ++i) {
3201        sum += blossomValue(i) * (blossomSize(i) / 2);
3202      }
3203      return sum;
3204    }
3205
3206    /// \brief Return the dual value (potential) of the given node.
3207    ///
3208    /// This function returns the dual value (potential) of the given node.
3209    ///
3210    /// \pre Either run() or start() must be called before using this function.
3211    Value nodeValue(const Node& n) const {
3212      return (*_node_potential)[n];
3213    }
3214
3215    /// \brief Return the number of the blossoms in the basis.
3216    ///
3217    /// This function returns the number of the blossoms in the basis.
3218    ///
3219    /// \pre Either run() or start() must be called before using this function.
3220    /// \see BlossomIt
3221    int blossomNum() const {
3222      return _blossom_potential.size();
3223    }
3224
3225    /// \brief Return the number of the nodes in the given blossom.
3226    ///
3227    /// This function returns the number of the nodes in the given blossom.
3228    ///
3229    /// \pre Either run() or start() must be called before using this function.
3230    /// \see BlossomIt
3231    int blossomSize(int k) const {
3232      return _blossom_potential[k].end - _blossom_potential[k].begin;
3233    }
3234
3235    /// \brief Return the dual value (ptential) of the given blossom.
3236    ///
3237    /// This function returns the dual value (ptential) of the given blossom.
3238    ///
3239    /// \pre Either run() or start() must be called before using this function.
3240    Value blossomValue(int k) const {
3241      return _blossom_potential[k].value;
3242    }
3243
3244    /// \brief Iterator for obtaining the nodes of a blossom.
3245    ///
3246    /// This class provides an iterator for obtaining the nodes of the
3247    /// given blossom. It lists a subset of the nodes.
3248    /// Before using this iterator, you must allocate a
3249    /// MaxWeightedPerfectMatching class and execute it.
3250    class BlossomIt {
3251    public:
3252
3253      /// \brief Constructor.
3254      ///
3255      /// Constructor to get the nodes of the given variable.
3256      ///
3257      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
3258      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
3259      /// must be called before initializing this iterator.
3260      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
3261        : _algorithm(&algorithm)
3262      {
3263        _index = _algorithm->_blossom_potential[variable].begin;
3264        _last = _algorithm->_blossom_potential[variable].end;
3265      }
3266
3267      /// \brief Conversion to \c Node.
3268      ///
3269      /// Conversion to \c Node.
3270      operator Node() const {
3271        return _algorithm->_blossom_node_list[_index];
3272      }
3273
3274      /// \brief Increment operator.
3275      ///
3276      /// Increment operator.
3277      BlossomIt& operator++() {
3278        ++_index;
3279        return *this;
3280      }
3281
3282      /// \brief Validity checking
3283      ///
3284      /// This function checks whether the iterator is invalid.
3285      bool operator==(Invalid) const { return _index == _last; }
3286
3287      /// \brief Validity checking
3288      ///
3289      /// This function checks whether the iterator is valid.
3290      bool operator!=(Invalid) const { return _index != _last; }
3291
3292    private:
3293      const MaxWeightedPerfectMatching* _algorithm;
3294      int _last;
3295      int _index;
3296    };
3297
3298    /// @}
3299
3300  };
3301
3302} //END OF NAMESPACE LEMON
3303
3304#endif //LEMON_MAX_MATCHING_H
Note: See TracBrowser for help on using the repository browser.