COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/max_matching.h @ 327:91d63b8b1a4c

Last change on this file since 327:91d63b8b1a4c was 327:91d63b8b1a4c, checked in by Balazs Dezso <deba@…>, 15 years ago

Several improvements in maximum matching algorithms

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