COIN-OR::LEMON - Graph Library

source: lemon/lemon/max_matching.h @ 613:e7017ec2d5cd

Last change on this file since 613:e7017ec2d5cd was 606:c5fd2d996909, checked in by Peter Kovacs <kpeter@…>, 15 years ago

Various doc improvements (#248)

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