COIN-OR::LEMON - Graph Library

source: lemon/lemon/max_matching.h @ 637:b61354458b59

Last change on this file since 637:b61354458b59 was 637:b61354458b59, checked in by Peter Kovacs <kpeter@…>, 11 years ago

Imporvements for the matching algorithms (#264)

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