COIN-OR::LEMON - Graph Library

source: lemon/lemon/max_matching.h @ 451:09e416d35896

Last change on this file since 451:09e416d35896 was 342:5ba887b7def4, checked in by Alpar Juttner <alpar@…>, 15 years ago

Doc improvements in lemon/max_matching.h

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