COIN-OR::LEMON - Graph Library

source: lemon/lemon/matching.h @ 948:636dadefe1e6

Last change on this file since 948:636dadefe1e6 was 947:0513ccfea967, checked in by Balazs Dezso <deba@…>, 11 years ago

General improvements in weighted matching algorithms (#314)

  • Fix include guard
  • Uniform handling of MATCHED and UNMATCHED blossoms
  • Prefer operations which decrease the number of trees
  • Fix improper use of '/='

The solved problems did not cause wrong solution.

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