COIN-OR::LEMON - Graph Library

source: lemon/lemon/matching.h @ 1337:4add05447ca0

Last change on this file since 1337:4add05447ca0 was 1270:dceba191c00d, checked in by Alpar Juttner <alpar@…>, 11 years ago

Apply unify-sources.sh to the source tree

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