COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/matching.h @ 871:86613aa28a0c

Last change on this file since 871:86613aa28a0c was 870:61120524af27, checked in by Balazs Dezso <deba@…>, 10 years ago

Fractional matching initialization of weighted matchings (#314)

File size: 108.1 KB
Line 
1/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library.
4 *
5 * Copyright (C) 2003-2009
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_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      if (!_node_potential) {
814        _node_potential = new NodePotential(_graph);
815      }
816      if (!_blossom_set) {
817        _blossom_index = new IntNodeMap(_graph);
818        _blossom_set = new BlossomSet(*_blossom_index);
819        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
820      }
821
822      if (!_node_index) {
823        _node_index = new IntNodeMap(_graph);
824        _node_heap_index = new IntArcMap(_graph);
825        _node_data = new RangeMap<NodeData>(_node_num,
826                                              NodeData(*_node_heap_index));
827      }
828
829      if (!_tree_set) {
830        _tree_set_index = new IntIntMap(_blossom_num);
831        _tree_set = new TreeSet(*_tree_set_index);
832      }
833      if (!_delta1) {
834        _delta1_index = new IntNodeMap(_graph);
835        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
836      }
837      if (!_delta2) {
838        _delta2_index = new IntIntMap(_blossom_num);
839        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
840      }
841      if (!_delta3) {
842        _delta3_index = new IntEdgeMap(_graph);
843        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
844      }
845      if (!_delta4) {
846        _delta4_index = new IntIntMap(_blossom_num);
847        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
848      }
849    }
850
851    void destroyStructures() {
852      if (_matching) {
853        delete _matching;
854      }
855      if (_node_potential) {
856        delete _node_potential;
857      }
858      if (_blossom_set) {
859        delete _blossom_index;
860        delete _blossom_set;
861        delete _blossom_data;
862      }
863
864      if (_node_index) {
865        delete _node_index;
866        delete _node_heap_index;
867        delete _node_data;
868      }
869
870      if (_tree_set) {
871        delete _tree_set_index;
872        delete _tree_set;
873      }
874      if (_delta1) {
875        delete _delta1_index;
876        delete _delta1;
877      }
878      if (_delta2) {
879        delete _delta2_index;
880        delete _delta2;
881      }
882      if (_delta3) {
883        delete _delta3_index;
884        delete _delta3;
885      }
886      if (_delta4) {
887        delete _delta4_index;
888        delete _delta4;
889      }
890    }
891
892    void matchedToEven(int blossom, int tree) {
893      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
894        _delta2->erase(blossom);
895      }
896
897      if (!_blossom_set->trivial(blossom)) {
898        (*_blossom_data)[blossom].pot -=
899          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
900      }
901
902      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
903           n != INVALID; ++n) {
904
905        _blossom_set->increase(n, std::numeric_limits<Value>::max());
906        int ni = (*_node_index)[n];
907
908        (*_node_data)[ni].heap.clear();
909        (*_node_data)[ni].heap_index.clear();
910
911        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
912
913        _delta1->push(n, (*_node_data)[ni].pot);
914
915        for (InArcIt e(_graph, n); e != INVALID; ++e) {
916          Node v = _graph.source(e);
917          int vb = _blossom_set->find(v);
918          int vi = (*_node_index)[v];
919
920          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
921            dualScale * _weight[e];
922
923          if ((*_blossom_data)[vb].status == EVEN) {
924            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
925              _delta3->push(e, rw / 2);
926            }
927          } else {
928            typename std::map<int, Arc>::iterator it =
929              (*_node_data)[vi].heap_index.find(tree);
930
931            if (it != (*_node_data)[vi].heap_index.end()) {
932              if ((*_node_data)[vi].heap[it->second] > rw) {
933                (*_node_data)[vi].heap.replace(it->second, e);
934                (*_node_data)[vi].heap.decrease(e, rw);
935                it->second = e;
936              }
937            } else {
938              (*_node_data)[vi].heap.push(e, rw);
939              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
940            }
941
942            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
943              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
944
945              if ((*_blossom_data)[vb].status == MATCHED) {
946                if (_delta2->state(vb) != _delta2->IN_HEAP) {
947                  _delta2->push(vb, _blossom_set->classPrio(vb) -
948                               (*_blossom_data)[vb].offset);
949                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
950                           (*_blossom_data)[vb].offset) {
951                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
952                                   (*_blossom_data)[vb].offset);
953                }
954              }
955            }
956          }
957        }
958      }
959      (*_blossom_data)[blossom].offset = 0;
960    }
961
962    void matchedToOdd(int blossom) {
963      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
964        _delta2->erase(blossom);
965      }
966      (*_blossom_data)[blossom].offset += _delta_sum;
967      if (!_blossom_set->trivial(blossom)) {
968        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
969                      (*_blossom_data)[blossom].offset);
970      }
971    }
972
973    void evenToMatched(int blossom, int tree) {
974      if (!_blossom_set->trivial(blossom)) {
975        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
976      }
977
978      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
979           n != INVALID; ++n) {
980        int ni = (*_node_index)[n];
981        (*_node_data)[ni].pot -= _delta_sum;
982
983        _delta1->erase(n);
984
985        for (InArcIt e(_graph, n); e != INVALID; ++e) {
986          Node v = _graph.source(e);
987          int vb = _blossom_set->find(v);
988          int vi = (*_node_index)[v];
989
990          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
991            dualScale * _weight[e];
992
993          if (vb == blossom) {
994            if (_delta3->state(e) == _delta3->IN_HEAP) {
995              _delta3->erase(e);
996            }
997          } else if ((*_blossom_data)[vb].status == EVEN) {
998
999            if (_delta3->state(e) == _delta3->IN_HEAP) {
1000              _delta3->erase(e);
1001            }
1002
1003            int vt = _tree_set->find(vb);
1004
1005            if (vt != tree) {
1006
1007              Arc r = _graph.oppositeArc(e);
1008
1009              typename std::map<int, Arc>::iterator it =
1010                (*_node_data)[ni].heap_index.find(vt);
1011
1012              if (it != (*_node_data)[ni].heap_index.end()) {
1013                if ((*_node_data)[ni].heap[it->second] > rw) {
1014                  (*_node_data)[ni].heap.replace(it->second, r);
1015                  (*_node_data)[ni].heap.decrease(r, rw);
1016                  it->second = r;
1017                }
1018              } else {
1019                (*_node_data)[ni].heap.push(r, rw);
1020                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1021              }
1022
1023              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1024                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1025
1026                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1027                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1028                               (*_blossom_data)[blossom].offset);
1029                } else if ((*_delta2)[blossom] >
1030                           _blossom_set->classPrio(blossom) -
1031                           (*_blossom_data)[blossom].offset){
1032                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1033                                   (*_blossom_data)[blossom].offset);
1034                }
1035              }
1036            }
1037          } else {
1038
1039            typename std::map<int, Arc>::iterator it =
1040              (*_node_data)[vi].heap_index.find(tree);
1041
1042            if (it != (*_node_data)[vi].heap_index.end()) {
1043              (*_node_data)[vi].heap.erase(it->second);
1044              (*_node_data)[vi].heap_index.erase(it);
1045              if ((*_node_data)[vi].heap.empty()) {
1046                _blossom_set->increase(v, std::numeric_limits<Value>::max());
1047              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1048                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1049              }
1050
1051              if ((*_blossom_data)[vb].status == MATCHED) {
1052                if (_blossom_set->classPrio(vb) ==
1053                    std::numeric_limits<Value>::max()) {
1054                  _delta2->erase(vb);
1055                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1056                           (*_blossom_data)[vb].offset) {
1057                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
1058                                   (*_blossom_data)[vb].offset);
1059                }
1060              }
1061            }
1062          }
1063        }
1064      }
1065    }
1066
1067    void oddToMatched(int blossom) {
1068      (*_blossom_data)[blossom].offset -= _delta_sum;
1069
1070      if (_blossom_set->classPrio(blossom) !=
1071          std::numeric_limits<Value>::max()) {
1072        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1073                      (*_blossom_data)[blossom].offset);
1074      }
1075
1076      if (!_blossom_set->trivial(blossom)) {
1077        _delta4->erase(blossom);
1078      }
1079    }
1080
1081    void oddToEven(int blossom, int tree) {
1082      if (!_blossom_set->trivial(blossom)) {
1083        _delta4->erase(blossom);
1084        (*_blossom_data)[blossom].pot -=
1085          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1086      }
1087
1088      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1089           n != INVALID; ++n) {
1090        int ni = (*_node_index)[n];
1091
1092        _blossom_set->increase(n, std::numeric_limits<Value>::max());
1093
1094        (*_node_data)[ni].heap.clear();
1095        (*_node_data)[ni].heap_index.clear();
1096        (*_node_data)[ni].pot +=
1097          2 * _delta_sum - (*_blossom_data)[blossom].offset;
1098
1099        _delta1->push(n, (*_node_data)[ni].pot);
1100
1101        for (InArcIt e(_graph, n); e != INVALID; ++e) {
1102          Node v = _graph.source(e);
1103          int vb = _blossom_set->find(v);
1104          int vi = (*_node_index)[v];
1105
1106          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1107            dualScale * _weight[e];
1108
1109          if ((*_blossom_data)[vb].status == EVEN) {
1110            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1111              _delta3->push(e, rw / 2);
1112            }
1113          } else {
1114
1115            typename std::map<int, Arc>::iterator it =
1116              (*_node_data)[vi].heap_index.find(tree);
1117
1118            if (it != (*_node_data)[vi].heap_index.end()) {
1119              if ((*_node_data)[vi].heap[it->second] > rw) {
1120                (*_node_data)[vi].heap.replace(it->second, e);
1121                (*_node_data)[vi].heap.decrease(e, rw);
1122                it->second = e;
1123              }
1124            } else {
1125              (*_node_data)[vi].heap.push(e, rw);
1126              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1127            }
1128
1129            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1130              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1131
1132              if ((*_blossom_data)[vb].status == MATCHED) {
1133                if (_delta2->state(vb) != _delta2->IN_HEAP) {
1134                  _delta2->push(vb, _blossom_set->classPrio(vb) -
1135                               (*_blossom_data)[vb].offset);
1136                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1137                           (*_blossom_data)[vb].offset) {
1138                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1139                                   (*_blossom_data)[vb].offset);
1140                }
1141              }
1142            }
1143          }
1144        }
1145      }
1146      (*_blossom_data)[blossom].offset = 0;
1147    }
1148
1149    void alternatePath(int even, int tree) {
1150      int odd;
1151
1152      evenToMatched(even, tree);
1153      (*_blossom_data)[even].status = MATCHED;
1154
1155      while ((*_blossom_data)[even].pred != INVALID) {
1156        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
1157        (*_blossom_data)[odd].status = MATCHED;
1158        oddToMatched(odd);
1159        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1160
1161        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
1162        (*_blossom_data)[even].status = MATCHED;
1163        evenToMatched(even, tree);
1164        (*_blossom_data)[even].next =
1165          _graph.oppositeArc((*_blossom_data)[odd].pred);
1166      }
1167
1168    }
1169
1170    void destroyTree(int tree) {
1171      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1172        if ((*_blossom_data)[b].status == EVEN) {
1173          (*_blossom_data)[b].status = MATCHED;
1174          evenToMatched(b, tree);
1175        } else if ((*_blossom_data)[b].status == ODD) {
1176          (*_blossom_data)[b].status = MATCHED;
1177          oddToMatched(b);
1178        }
1179      }
1180      _tree_set->eraseClass(tree);
1181    }
1182
1183
1184    void unmatchNode(const Node& node) {
1185      int blossom = _blossom_set->find(node);
1186      int tree = _tree_set->find(blossom);
1187
1188      alternatePath(blossom, tree);
1189      destroyTree(tree);
1190
1191      (*_blossom_data)[blossom].base = node;
1192      (*_blossom_data)[blossom].next = INVALID;
1193    }
1194
1195    void augmentOnEdge(const Edge& edge) {
1196
1197      int left = _blossom_set->find(_graph.u(edge));
1198      int right = _blossom_set->find(_graph.v(edge));
1199
1200      int left_tree = _tree_set->find(left);
1201      alternatePath(left, left_tree);
1202      destroyTree(left_tree);
1203
1204      int right_tree = _tree_set->find(right);
1205      alternatePath(right, right_tree);
1206      destroyTree(right_tree);
1207
1208      (*_blossom_data)[left].next = _graph.direct(edge, true);
1209      (*_blossom_data)[right].next = _graph.direct(edge, false);
1210    }
1211
1212    void augmentOnArc(const Arc& arc) {
1213
1214      int left = _blossom_set->find(_graph.source(arc));
1215      int right = _blossom_set->find(_graph.target(arc));
1216
1217      (*_blossom_data)[left].status = MATCHED;
1218
1219      int right_tree = _tree_set->find(right);
1220      alternatePath(right, right_tree);
1221      destroyTree(right_tree);
1222
1223      (*_blossom_data)[left].next = arc;
1224      (*_blossom_data)[right].next = _graph.oppositeArc(arc);
1225    }
1226
1227    void extendOnArc(const Arc& arc) {
1228      int base = _blossom_set->find(_graph.target(arc));
1229      int tree = _tree_set->find(base);
1230
1231      int odd = _blossom_set->find(_graph.source(arc));
1232      _tree_set->insert(odd, tree);
1233      (*_blossom_data)[odd].status = ODD;
1234      matchedToOdd(odd);
1235      (*_blossom_data)[odd].pred = arc;
1236
1237      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
1238      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1239      _tree_set->insert(even, tree);
1240      (*_blossom_data)[even].status = EVEN;
1241      matchedToEven(even, tree);
1242    }
1243
1244    void shrinkOnEdge(const Edge& edge, int tree) {
1245      int nca = -1;
1246      std::vector<int> left_path, right_path;
1247
1248      {
1249        std::set<int> left_set, right_set;
1250        int left = _blossom_set->find(_graph.u(edge));
1251        left_path.push_back(left);
1252        left_set.insert(left);
1253
1254        int right = _blossom_set->find(_graph.v(edge));
1255        right_path.push_back(right);
1256        right_set.insert(right);
1257
1258        while (true) {
1259
1260          if ((*_blossom_data)[left].pred == INVALID) break;
1261
1262          left =
1263            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1264          left_path.push_back(left);
1265          left =
1266            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1267          left_path.push_back(left);
1268
1269          left_set.insert(left);
1270
1271          if (right_set.find(left) != right_set.end()) {
1272            nca = left;
1273            break;
1274          }
1275
1276          if ((*_blossom_data)[right].pred == INVALID) break;
1277
1278          right =
1279            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1280          right_path.push_back(right);
1281          right =
1282            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1283          right_path.push_back(right);
1284
1285          right_set.insert(right);
1286
1287          if (left_set.find(right) != left_set.end()) {
1288            nca = right;
1289            break;
1290          }
1291
1292        }
1293
1294        if (nca == -1) {
1295          if ((*_blossom_data)[left].pred == INVALID) {
1296            nca = right;
1297            while (left_set.find(nca) == left_set.end()) {
1298              nca =
1299                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1300              right_path.push_back(nca);
1301              nca =
1302                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1303              right_path.push_back(nca);
1304            }
1305          } else {
1306            nca = left;
1307            while (right_set.find(nca) == right_set.end()) {
1308              nca =
1309                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1310              left_path.push_back(nca);
1311              nca =
1312                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1313              left_path.push_back(nca);
1314            }
1315          }
1316        }
1317      }
1318
1319      std::vector<int> subblossoms;
1320      Arc prev;
1321
1322      prev = _graph.direct(edge, true);
1323      for (int i = 0; left_path[i] != nca; i += 2) {
1324        subblossoms.push_back(left_path[i]);
1325        (*_blossom_data)[left_path[i]].next = prev;
1326        _tree_set->erase(left_path[i]);
1327
1328        subblossoms.push_back(left_path[i + 1]);
1329        (*_blossom_data)[left_path[i + 1]].status = EVEN;
1330        oddToEven(left_path[i + 1], tree);
1331        _tree_set->erase(left_path[i + 1]);
1332        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
1333      }
1334
1335      int k = 0;
1336      while (right_path[k] != nca) ++k;
1337
1338      subblossoms.push_back(nca);
1339      (*_blossom_data)[nca].next = prev;
1340
1341      for (int i = k - 2; i >= 0; i -= 2) {
1342        subblossoms.push_back(right_path[i + 1]);
1343        (*_blossom_data)[right_path[i + 1]].status = EVEN;
1344        oddToEven(right_path[i + 1], tree);
1345        _tree_set->erase(right_path[i + 1]);
1346
1347        (*_blossom_data)[right_path[i + 1]].next =
1348          (*_blossom_data)[right_path[i + 1]].pred;
1349
1350        subblossoms.push_back(right_path[i]);
1351        _tree_set->erase(right_path[i]);
1352      }
1353
1354      int surface =
1355        _blossom_set->join(subblossoms.begin(), subblossoms.end());
1356
1357      for (int i = 0; i < int(subblossoms.size()); ++i) {
1358        if (!_blossom_set->trivial(subblossoms[i])) {
1359          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1360        }
1361        (*_blossom_data)[subblossoms[i]].status = MATCHED;
1362      }
1363
1364      (*_blossom_data)[surface].pot = -2 * _delta_sum;
1365      (*_blossom_data)[surface].offset = 0;
1366      (*_blossom_data)[surface].status = EVEN;
1367      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1368      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1369
1370      _tree_set->insert(surface, tree);
1371      _tree_set->erase(nca);
1372    }
1373
1374    void splitBlossom(int blossom) {
1375      Arc next = (*_blossom_data)[blossom].next;
1376      Arc pred = (*_blossom_data)[blossom].pred;
1377
1378      int tree = _tree_set->find(blossom);
1379
1380      (*_blossom_data)[blossom].status = MATCHED;
1381      oddToMatched(blossom);
1382      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1383        _delta2->erase(blossom);
1384      }
1385
1386      std::vector<int> subblossoms;
1387      _blossom_set->split(blossom, std::back_inserter(subblossoms));
1388
1389      Value offset = (*_blossom_data)[blossom].offset;
1390      int b = _blossom_set->find(_graph.source(pred));
1391      int d = _blossom_set->find(_graph.source(next));
1392
1393      int ib = -1, id = -1;
1394      for (int i = 0; i < int(subblossoms.size()); ++i) {
1395        if (subblossoms[i] == b) ib = i;
1396        if (subblossoms[i] == d) id = i;
1397
1398        (*_blossom_data)[subblossoms[i]].offset = offset;
1399        if (!_blossom_set->trivial(subblossoms[i])) {
1400          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1401        }
1402        if (_blossom_set->classPrio(subblossoms[i]) !=
1403            std::numeric_limits<Value>::max()) {
1404          _delta2->push(subblossoms[i],
1405                        _blossom_set->classPrio(subblossoms[i]) -
1406                        (*_blossom_data)[subblossoms[i]].offset);
1407        }
1408      }
1409
1410      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1411        for (int i = (id + 1) % subblossoms.size();
1412             i != ib; i = (i + 2) % subblossoms.size()) {
1413          int sb = subblossoms[i];
1414          int tb = subblossoms[(i + 1) % subblossoms.size()];
1415          (*_blossom_data)[sb].next =
1416            _graph.oppositeArc((*_blossom_data)[tb].next);
1417        }
1418
1419        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1420          int sb = subblossoms[i];
1421          int tb = subblossoms[(i + 1) % subblossoms.size()];
1422          int ub = subblossoms[(i + 2) % subblossoms.size()];
1423
1424          (*_blossom_data)[sb].status = ODD;
1425          matchedToOdd(sb);
1426          _tree_set->insert(sb, tree);
1427          (*_blossom_data)[sb].pred = pred;
1428          (*_blossom_data)[sb].next =
1429            _graph.oppositeArc((*_blossom_data)[tb].next);
1430
1431          pred = (*_blossom_data)[ub].next;
1432
1433          (*_blossom_data)[tb].status = EVEN;
1434          matchedToEven(tb, tree);
1435          _tree_set->insert(tb, tree);
1436          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1437        }
1438
1439        (*_blossom_data)[subblossoms[id]].status = ODD;
1440        matchedToOdd(subblossoms[id]);
1441        _tree_set->insert(subblossoms[id], tree);
1442        (*_blossom_data)[subblossoms[id]].next = next;
1443        (*_blossom_data)[subblossoms[id]].pred = pred;
1444
1445      } else {
1446
1447        for (int i = (ib + 1) % subblossoms.size();
1448             i != id; i = (i + 2) % subblossoms.size()) {
1449          int sb = subblossoms[i];
1450          int tb = subblossoms[(i + 1) % subblossoms.size()];
1451          (*_blossom_data)[sb].next =
1452            _graph.oppositeArc((*_blossom_data)[tb].next);
1453        }
1454
1455        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1456          int sb = subblossoms[i];
1457          int tb = subblossoms[(i + 1) % subblossoms.size()];
1458          int ub = subblossoms[(i + 2) % subblossoms.size()];
1459
1460          (*_blossom_data)[sb].status = ODD;
1461          matchedToOdd(sb);
1462          _tree_set->insert(sb, tree);
1463          (*_blossom_data)[sb].next = next;
1464          (*_blossom_data)[sb].pred =
1465            _graph.oppositeArc((*_blossom_data)[tb].next);
1466
1467          (*_blossom_data)[tb].status = EVEN;
1468          matchedToEven(tb, tree);
1469          _tree_set->insert(tb, tree);
1470          (*_blossom_data)[tb].pred =
1471            (*_blossom_data)[tb].next =
1472            _graph.oppositeArc((*_blossom_data)[ub].next);
1473          next = (*_blossom_data)[ub].next;
1474        }
1475
1476        (*_blossom_data)[subblossoms[ib]].status = ODD;
1477        matchedToOdd(subblossoms[ib]);
1478        _tree_set->insert(subblossoms[ib], tree);
1479        (*_blossom_data)[subblossoms[ib]].next = next;
1480        (*_blossom_data)[subblossoms[ib]].pred = pred;
1481      }
1482      _tree_set->erase(blossom);
1483    }
1484
1485    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1486      if (_blossom_set->trivial(blossom)) {
1487        int bi = (*_node_index)[base];
1488        Value pot = (*_node_data)[bi].pot;
1489
1490        (*_matching)[base] = matching;
1491        _blossom_node_list.push_back(base);
1492        (*_node_potential)[base] = pot;
1493      } else {
1494
1495        Value pot = (*_blossom_data)[blossom].pot;
1496        int bn = _blossom_node_list.size();
1497
1498        std::vector<int> subblossoms;
1499        _blossom_set->split(blossom, std::back_inserter(subblossoms));
1500        int b = _blossom_set->find(base);
1501        int ib = -1;
1502        for (int i = 0; i < int(subblossoms.size()); ++i) {
1503          if (subblossoms[i] == b) { ib = i; break; }
1504        }
1505
1506        for (int i = 1; i < int(subblossoms.size()); i += 2) {
1507          int sb = subblossoms[(ib + i) % subblossoms.size()];
1508          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1509
1510          Arc m = (*_blossom_data)[tb].next;
1511          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1512          extractBlossom(tb, _graph.source(m), m);
1513        }
1514        extractBlossom(subblossoms[ib], base, matching);
1515
1516        int en = _blossom_node_list.size();
1517
1518        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1519      }
1520    }
1521
1522    void extractMatching() {
1523      std::vector<int> blossoms;
1524      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1525        blossoms.push_back(c);
1526      }
1527
1528      for (int i = 0; i < int(blossoms.size()); ++i) {
1529        if ((*_blossom_data)[blossoms[i]].next != INVALID) {
1530
1531          Value offset = (*_blossom_data)[blossoms[i]].offset;
1532          (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1533          for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1534               n != INVALID; ++n) {
1535            (*_node_data)[(*_node_index)[n]].pot -= offset;
1536          }
1537
1538          Arc matching = (*_blossom_data)[blossoms[i]].next;
1539          Node base = _graph.source(matching);
1540          extractBlossom(blossoms[i], base, matching);
1541        } else {
1542          Node base = (*_blossom_data)[blossoms[i]].base;
1543          extractBlossom(blossoms[i], base, INVALID);
1544        }
1545      }
1546    }
1547
1548  public:
1549
1550    /// \brief Constructor
1551    ///
1552    /// Constructor.
1553    MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1554      : _graph(graph), _weight(weight), _matching(0),
1555        _node_potential(0), _blossom_potential(), _blossom_node_list(),
1556        _node_num(0), _blossom_num(0),
1557
1558        _blossom_index(0), _blossom_set(0), _blossom_data(0),
1559        _node_index(0), _node_heap_index(0), _node_data(0),
1560        _tree_set_index(0), _tree_set(0),
1561
1562        _delta1_index(0), _delta1(0),
1563        _delta2_index(0), _delta2(0),
1564        _delta3_index(0), _delta3(0),
1565        _delta4_index(0), _delta4(0),
1566
1567        _delta_sum(), _unmatched(0),
1568
1569        _fractional(0)
1570    {}
1571
1572    ~MaxWeightedMatching() {
1573      destroyStructures();
1574      if (_fractional) {
1575        delete _fractional;
1576      }
1577    }
1578
1579    /// \name Execution Control
1580    /// The simplest way to execute the algorithm is to use the
1581    /// \ref run() member function.
1582
1583    ///@{
1584
1585    /// \brief Initialize the algorithm
1586    ///
1587    /// This function initializes the algorithm.
1588    void init() {
1589      createStructures();
1590
1591      for (ArcIt e(_graph); e != INVALID; ++e) {
1592        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1593      }
1594      for (NodeIt n(_graph); n != INVALID; ++n) {
1595        (*_delta1_index)[n] = _delta1->PRE_HEAP;
1596      }
1597      for (EdgeIt e(_graph); e != INVALID; ++e) {
1598        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1599      }
1600      for (int i = 0; i < _blossom_num; ++i) {
1601        (*_delta2_index)[i] = _delta2->PRE_HEAP;
1602        (*_delta4_index)[i] = _delta4->PRE_HEAP;
1603      }
1604
1605      _unmatched = _node_num;
1606
1607      int index = 0;
1608      for (NodeIt n(_graph); n != INVALID; ++n) {
1609        Value max = 0;
1610        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1611          if (_graph.target(e) == n) continue;
1612          if ((dualScale * _weight[e]) / 2 > max) {
1613            max = (dualScale * _weight[e]) / 2;
1614          }
1615        }
1616        (*_node_index)[n] = index;
1617        (*_node_data)[index].pot = max;
1618        _delta1->push(n, max);
1619        int blossom =
1620          _blossom_set->insert(n, std::numeric_limits<Value>::max());
1621
1622        _tree_set->insert(blossom);
1623
1624        (*_blossom_data)[blossom].status = EVEN;
1625        (*_blossom_data)[blossom].pred = INVALID;
1626        (*_blossom_data)[blossom].next = INVALID;
1627        (*_blossom_data)[blossom].pot = 0;
1628        (*_blossom_data)[blossom].offset = 0;
1629        ++index;
1630      }
1631      for (EdgeIt e(_graph); e != INVALID; ++e) {
1632        int si = (*_node_index)[_graph.u(e)];
1633        int ti = (*_node_index)[_graph.v(e)];
1634        if (_graph.u(e) != _graph.v(e)) {
1635          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1636                            dualScale * _weight[e]) / 2);
1637        }
1638      }
1639    }
1640
1641    /// \brief Initialize the algorithm with fractional matching
1642    ///
1643    /// This function initializes the algorithm with a fractional
1644    /// matching. This initialization is also called jumpstart heuristic.
1645    void fractionalInit() {
1646      createStructures();
1647     
1648      if (_fractional == 0) {
1649        _fractional = new FractionalMatching(_graph, _weight, false);
1650      }
1651      _fractional->run();
1652
1653      for (ArcIt e(_graph); e != INVALID; ++e) {
1654        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1655      }
1656      for (NodeIt n(_graph); n != INVALID; ++n) {
1657        (*_delta1_index)[n] = _delta1->PRE_HEAP;
1658      }
1659      for (EdgeIt e(_graph); e != INVALID; ++e) {
1660        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1661      }
1662      for (int i = 0; i < _blossom_num; ++i) {
1663        (*_delta2_index)[i] = _delta2->PRE_HEAP;
1664        (*_delta4_index)[i] = _delta4->PRE_HEAP;
1665      }
1666
1667      _unmatched = 0;
1668
1669      int index = 0;
1670      for (NodeIt n(_graph); n != INVALID; ++n) {
1671        Value pot = _fractional->nodeValue(n);
1672        (*_node_index)[n] = index;
1673        (*_node_data)[index].pot = pot;
1674        int blossom =
1675          _blossom_set->insert(n, std::numeric_limits<Value>::max());
1676
1677        (*_blossom_data)[blossom].status = MATCHED;
1678        (*_blossom_data)[blossom].pred = INVALID;
1679        (*_blossom_data)[blossom].next = _fractional->matching(n);
1680        if (_fractional->matching(n) == INVALID) {
1681          (*_blossom_data)[blossom].base = n;
1682        }
1683        (*_blossom_data)[blossom].pot = 0;
1684        (*_blossom_data)[blossom].offset = 0;
1685        ++index;
1686      }
1687
1688      typename Graph::template NodeMap<bool> processed(_graph, false);
1689      for (NodeIt n(_graph); n != INVALID; ++n) {
1690        if (processed[n]) continue;
1691        processed[n] = true;
1692        if (_fractional->matching(n) == INVALID) continue;
1693        int num = 1;
1694        Node v = _graph.target(_fractional->matching(n));
1695        while (n != v) {
1696          processed[v] = true;
1697          v = _graph.target(_fractional->matching(v));
1698          ++num;
1699        }
1700
1701        if (num % 2 == 1) {
1702          std::vector<int> subblossoms(num);
1703
1704          subblossoms[--num] = _blossom_set->find(n);
1705          _delta1->push(n, _fractional->nodeValue(n));
1706          v = _graph.target(_fractional->matching(n));
1707          while (n != v) {
1708            subblossoms[--num] = _blossom_set->find(v);
1709            _delta1->push(v, _fractional->nodeValue(v));
1710            v = _graph.target(_fractional->matching(v));           
1711          }
1712         
1713          int surface =
1714            _blossom_set->join(subblossoms.begin(), subblossoms.end());
1715          (*_blossom_data)[surface].status = EVEN;
1716          (*_blossom_data)[surface].pred = INVALID;
1717          (*_blossom_data)[surface].next = INVALID;
1718          (*_blossom_data)[surface].pot = 0;
1719          (*_blossom_data)[surface].offset = 0;
1720         
1721          _tree_set->insert(surface);
1722          ++_unmatched;
1723        }
1724      }
1725
1726      for (EdgeIt e(_graph); e != INVALID; ++e) {
1727        int si = (*_node_index)[_graph.u(e)];
1728        int sb = _blossom_set->find(_graph.u(e));
1729        int ti = (*_node_index)[_graph.v(e)];
1730        int tb = _blossom_set->find(_graph.v(e));
1731        if ((*_blossom_data)[sb].status == EVEN &&
1732            (*_blossom_data)[tb].status == EVEN && sb != tb) {
1733          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1734                            dualScale * _weight[e]) / 2);
1735        }
1736      }
1737
1738      for (NodeIt n(_graph); n != INVALID; ++n) {
1739        int nb = _blossom_set->find(n);
1740        if ((*_blossom_data)[nb].status != MATCHED) continue;
1741        int ni = (*_node_index)[n];
1742
1743        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1744          Node v = _graph.target(e);
1745          int vb = _blossom_set->find(v);
1746          int vi = (*_node_index)[v];
1747
1748          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1749            dualScale * _weight[e];
1750
1751          if ((*_blossom_data)[vb].status == EVEN) {
1752
1753            int vt = _tree_set->find(vb);
1754
1755            typename std::map<int, Arc>::iterator it =
1756              (*_node_data)[ni].heap_index.find(vt);
1757
1758            if (it != (*_node_data)[ni].heap_index.end()) {
1759              if ((*_node_data)[ni].heap[it->second] > rw) {
1760                (*_node_data)[ni].heap.replace(it->second, e);
1761                (*_node_data)[ni].heap.decrease(e, rw);
1762                it->second = e;
1763              }
1764            } else {
1765              (*_node_data)[ni].heap.push(e, rw);
1766              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
1767            }
1768          }
1769        }
1770           
1771        if (!(*_node_data)[ni].heap.empty()) {
1772          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1773          _delta2->push(nb, _blossom_set->classPrio(nb));
1774        }
1775      }
1776    }
1777
1778    /// \brief Start the algorithm
1779    ///
1780    /// This function starts the algorithm.
1781    ///
1782    /// \pre \ref init() or \ref fractionalInit() must be called
1783    /// before using this function.
1784    void start() {
1785      enum OpType {
1786        D1, D2, D3, D4
1787      };
1788
1789      while (_unmatched > 0) {
1790        Value d1 = !_delta1->empty() ?
1791          _delta1->prio() : std::numeric_limits<Value>::max();
1792
1793        Value d2 = !_delta2->empty() ?
1794          _delta2->prio() : std::numeric_limits<Value>::max();
1795
1796        Value d3 = !_delta3->empty() ?
1797          _delta3->prio() : std::numeric_limits<Value>::max();
1798
1799        Value d4 = !_delta4->empty() ?
1800          _delta4->prio() : std::numeric_limits<Value>::max();
1801
1802        _delta_sum = d3; OpType ot = D3;
1803        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1804        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1805        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1806
1807        switch (ot) {
1808        case D1:
1809          {
1810            Node n = _delta1->top();
1811            unmatchNode(n);
1812            --_unmatched;
1813          }
1814          break;
1815        case D2:
1816          {
1817            int blossom = _delta2->top();
1818            Node n = _blossom_set->classTop(blossom);
1819            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
1820            if ((*_blossom_data)[blossom].next == INVALID) {
1821              augmentOnArc(a);
1822              --_unmatched;
1823            } else {
1824              extendOnArc(a);
1825            }
1826          }
1827          break;
1828        case D3:
1829          {
1830            Edge e = _delta3->top();
1831
1832            int left_blossom = _blossom_set->find(_graph.u(e));
1833            int right_blossom = _blossom_set->find(_graph.v(e));
1834
1835            if (left_blossom == right_blossom) {
1836              _delta3->pop();
1837            } else {
1838              int left_tree = _tree_set->find(left_blossom);
1839              int right_tree = _tree_set->find(right_blossom);
1840
1841              if (left_tree == right_tree) {
1842                shrinkOnEdge(e, left_tree);
1843              } else {
1844                augmentOnEdge(e);
1845                _unmatched -= 2;
1846              }
1847            }
1848          } break;
1849        case D4:
1850          splitBlossom(_delta4->top());
1851          break;
1852        }
1853      }
1854      extractMatching();
1855    }
1856
1857    /// \brief Run the algorithm.
1858    ///
1859    /// This method runs the \c %MaxWeightedMatching algorithm.
1860    ///
1861    /// \note mwm.run() is just a shortcut of the following code.
1862    /// \code
1863    ///   mwm.fractionalInit();
1864    ///   mwm.start();
1865    /// \endcode
1866    void run() {
1867      fractionalInit();
1868      start();
1869    }
1870
1871    /// @}
1872
1873    /// \name Primal Solution
1874    /// Functions to get the primal solution, i.e. the maximum weighted
1875    /// matching.\n
1876    /// Either \ref run() or \ref start() function should be called before
1877    /// using them.
1878
1879    /// @{
1880
1881    /// \brief Return the weight of the matching.
1882    ///
1883    /// This function returns the weight of the found matching.
1884    ///
1885    /// \pre Either run() or start() must be called before using this function.
1886    Value matchingWeight() const {
1887      Value sum = 0;
1888      for (NodeIt n(_graph); n != INVALID; ++n) {
1889        if ((*_matching)[n] != INVALID) {
1890          sum += _weight[(*_matching)[n]];
1891        }
1892      }
1893      return sum / 2;
1894    }
1895
1896    /// \brief Return the size (cardinality) of the matching.
1897    ///
1898    /// This function returns the size (cardinality) of the found matching.
1899    ///
1900    /// \pre Either run() or start() must be called before using this function.
1901    int matchingSize() const {
1902      int num = 0;
1903      for (NodeIt n(_graph); n != INVALID; ++n) {
1904        if ((*_matching)[n] != INVALID) {
1905          ++num;
1906        }
1907      }
1908      return num /= 2;
1909    }
1910
1911    /// \brief Return \c true if the given edge is in the matching.
1912    ///
1913    /// This function returns \c true if the given edge is in the found
1914    /// matching.
1915    ///
1916    /// \pre Either run() or start() must be called before using this function.
1917    bool matching(const Edge& edge) const {
1918      return edge == (*_matching)[_graph.u(edge)];
1919    }
1920
1921    /// \brief Return the matching arc (or edge) incident to the given node.
1922    ///
1923    /// This function returns the matching arc (or edge) incident to the
1924    /// given node in the found matching or \c INVALID if the node is
1925    /// not covered by the matching.
1926    ///
1927    /// \pre Either run() or start() must be called before using this function.
1928    Arc matching(const Node& node) const {
1929      return (*_matching)[node];
1930    }
1931
1932    /// \brief Return a const reference to the matching map.
1933    ///
1934    /// This function returns a const reference to a node map that stores
1935    /// the matching arc (or edge) incident to each node.
1936    const MatchingMap& matchingMap() const {
1937      return *_matching;
1938    }
1939
1940    /// \brief Return the mate of the given node.
1941    ///
1942    /// This function returns the mate of the given node in the found
1943    /// matching or \c INVALID if the node is not covered by the matching.
1944    ///
1945    /// \pre Either run() or start() must be called before using this function.
1946    Node mate(const Node& node) const {
1947      return (*_matching)[node] != INVALID ?
1948        _graph.target((*_matching)[node]) : INVALID;
1949    }
1950
1951    /// @}
1952
1953    /// \name Dual Solution
1954    /// Functions to get the dual solution.\n
1955    /// Either \ref run() or \ref start() function should be called before
1956    /// using them.
1957
1958    /// @{
1959
1960    /// \brief Return the value of the dual solution.
1961    ///
1962    /// This function returns the value of the dual solution.
1963    /// It should be equal to the primal value scaled by \ref dualScale
1964    /// "dual scale".
1965    ///
1966    /// \pre Either run() or start() must be called before using this function.
1967    Value dualValue() const {
1968      Value sum = 0;
1969      for (NodeIt n(_graph); n != INVALID; ++n) {
1970        sum += nodeValue(n);
1971      }
1972      for (int i = 0; i < blossomNum(); ++i) {
1973        sum += blossomValue(i) * (blossomSize(i) / 2);
1974      }
1975      return sum;
1976    }
1977
1978    /// \brief Return the dual value (potential) of the given node.
1979    ///
1980    /// This function returns the dual value (potential) of the given node.
1981    ///
1982    /// \pre Either run() or start() must be called before using this function.
1983    Value nodeValue(const Node& n) const {
1984      return (*_node_potential)[n];
1985    }
1986
1987    /// \brief Return the number of the blossoms in the basis.
1988    ///
1989    /// This function returns the number of the blossoms in the basis.
1990    ///
1991    /// \pre Either run() or start() must be called before using this function.
1992    /// \see BlossomIt
1993    int blossomNum() const {
1994      return _blossom_potential.size();
1995    }
1996
1997    /// \brief Return the number of the nodes in the given blossom.
1998    ///
1999    /// This function returns the number of the nodes in the given blossom.
2000    ///
2001    /// \pre Either run() or start() must be called before using this function.
2002    /// \see BlossomIt
2003    int blossomSize(int k) const {
2004      return _blossom_potential[k].end - _blossom_potential[k].begin;
2005    }
2006
2007    /// \brief Return the dual value (ptential) of the given blossom.
2008    ///
2009    /// This function returns the dual value (ptential) of the given blossom.
2010    ///
2011    /// \pre Either run() or start() must be called before using this function.
2012    Value blossomValue(int k) const {
2013      return _blossom_potential[k].value;
2014    }
2015
2016    /// \brief Iterator for obtaining the nodes of a blossom.
2017    ///
2018    /// This class provides an iterator for obtaining the nodes of the
2019    /// given blossom. It lists a subset of the nodes.
2020    /// Before using this iterator, you must allocate a
2021    /// MaxWeightedMatching class and execute it.
2022    class BlossomIt {
2023    public:
2024
2025      /// \brief Constructor.
2026      ///
2027      /// Constructor to get the nodes of the given variable.
2028      ///
2029      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
2030      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
2031      /// called before initializing this iterator.
2032      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
2033        : _algorithm(&algorithm)
2034      {
2035        _index = _algorithm->_blossom_potential[variable].begin;
2036        _last = _algorithm->_blossom_potential[variable].end;
2037      }
2038
2039      /// \brief Conversion to \c Node.
2040      ///
2041      /// Conversion to \c Node.
2042      operator Node() const {
2043        return _algorithm->_blossom_node_list[_index];
2044      }
2045
2046      /// \brief Increment operator.
2047      ///
2048      /// Increment operator.
2049      BlossomIt& operator++() {
2050        ++_index;
2051        return *this;
2052      }
2053
2054      /// \brief Validity checking
2055      ///
2056      /// Checks whether the iterator is invalid.
2057      bool operator==(Invalid) const { return _index == _last; }
2058
2059      /// \brief Validity checking
2060      ///
2061      /// Checks whether the iterator is valid.
2062      bool operator!=(Invalid) const { return _index != _last; }
2063
2064    private:
2065      const MaxWeightedMatching* _algorithm;
2066      int _last;
2067      int _index;
2068    };
2069
2070    /// @}
2071
2072  };
2073
2074  /// \ingroup matching
2075  ///
2076  /// \brief Weighted perfect matching in general graphs
2077  ///
2078  /// This class provides an efficient implementation of Edmond's
2079  /// maximum weighted perfect matching algorithm. The implementation
2080  /// is based on extensive use of priority queues and provides
2081  /// \f$O(nm\log n)\f$ time complexity.
2082  ///
2083  /// The maximum weighted perfect matching problem is to find a subset of
2084  /// the edges in an undirected graph with maximum overall weight for which
2085  /// each node has exactly one incident edge.
2086  /// It can be formulated with the following linear program.
2087  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
2088  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
2089      \quad \forall B\in\mathcal{O}\f] */
2090  /// \f[x_e \ge 0\quad \forall e\in E\f]
2091  /// \f[\max \sum_{e\in E}x_ew_e\f]
2092  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
2093  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
2094  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
2095  /// subsets of the nodes.
2096  ///
2097  /// The algorithm calculates an optimal matching and a proof of the
2098  /// optimality. The solution of the dual problem can be used to check
2099  /// the result of the algorithm. The dual linear problem is the
2100  /// following.
2101  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
2102      w_{uv} \quad \forall uv\in E\f] */
2103  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
2104  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
2105      \frac{\vert B \vert - 1}{2}z_B\f] */
2106  ///
2107  /// The algorithm can be executed with the run() function.
2108  /// After it the matching (the primal solution) and the dual solution
2109  /// can be obtained using the query functions and the
2110  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
2111  /// which is able to iterate on the nodes of a blossom.
2112  /// If the value type is integer, then the dual solution is multiplied
2113  /// by \ref MaxWeightedMatching::dualScale "4".
2114  ///
2115  /// \tparam GR The undirected graph type the algorithm runs on.
2116  /// \tparam WM The type edge weight map. The default type is
2117  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
2118#ifdef DOXYGEN
2119  template <typename GR, typename WM>
2120#else
2121  template <typename GR,
2122            typename WM = typename GR::template EdgeMap<int> >
2123#endif
2124  class MaxWeightedPerfectMatching {
2125  public:
2126
2127    /// The graph type of the algorithm
2128    typedef GR Graph;
2129    /// The type of the edge weight map
2130    typedef WM WeightMap;
2131    /// The value type of the edge weights
2132    typedef typename WeightMap::Value Value;
2133
2134    /// \brief Scaling factor for dual solution
2135    ///
2136    /// Scaling factor for dual solution, it is equal to 4 or 1
2137    /// according to the value type.
2138    static const int dualScale =
2139      std::numeric_limits<Value>::is_integer ? 4 : 1;
2140
2141    /// The type of the matching map
2142    typedef typename Graph::template NodeMap<typename Graph::Arc>
2143    MatchingMap;
2144
2145  private:
2146
2147    TEMPLATE_GRAPH_TYPEDEFS(Graph);
2148
2149    typedef typename Graph::template NodeMap<Value> NodePotential;
2150    typedef std::vector<Node> BlossomNodeList;
2151
2152    struct BlossomVariable {
2153      int begin, end;
2154      Value value;
2155
2156      BlossomVariable(int _begin, int _end, Value _value)
2157        : begin(_begin), end(_end), value(_value) {}
2158
2159    };
2160
2161    typedef std::vector<BlossomVariable> BlossomPotential;
2162
2163    const Graph& _graph;
2164    const WeightMap& _weight;
2165
2166    MatchingMap* _matching;
2167
2168    NodePotential* _node_potential;
2169
2170    BlossomPotential _blossom_potential;
2171    BlossomNodeList _blossom_node_list;
2172
2173    int _node_num;
2174    int _blossom_num;
2175
2176    typedef RangeMap<int> IntIntMap;
2177
2178    enum Status {
2179      EVEN = -1, MATCHED = 0, ODD = 1
2180    };
2181
2182    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
2183    struct BlossomData {
2184      int tree;
2185      Status status;
2186      Arc pred, next;
2187      Value pot, offset;
2188    };
2189
2190    IntNodeMap *_blossom_index;
2191    BlossomSet *_blossom_set;
2192    RangeMap<BlossomData>* _blossom_data;
2193
2194    IntNodeMap *_node_index;
2195    IntArcMap *_node_heap_index;
2196
2197    struct NodeData {
2198
2199      NodeData(IntArcMap& node_heap_index)
2200        : heap(node_heap_index) {}
2201
2202      int blossom;
2203      Value pot;
2204      BinHeap<Value, IntArcMap> heap;
2205      std::map<int, Arc> heap_index;
2206
2207      int tree;
2208    };
2209
2210    RangeMap<NodeData>* _node_data;
2211
2212    typedef ExtendFindEnum<IntIntMap> TreeSet;
2213
2214    IntIntMap *_tree_set_index;
2215    TreeSet *_tree_set;
2216
2217    IntIntMap *_delta2_index;
2218    BinHeap<Value, IntIntMap> *_delta2;
2219
2220    IntEdgeMap *_delta3_index;
2221    BinHeap<Value, IntEdgeMap> *_delta3;
2222
2223    IntIntMap *_delta4_index;
2224    BinHeap<Value, IntIntMap> *_delta4;
2225
2226    Value _delta_sum;
2227    int _unmatched;
2228
2229    typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
2230    FractionalMatching;
2231    FractionalMatching *_fractional;
2232
2233    void createStructures() {
2234      _node_num = countNodes(_graph);
2235      _blossom_num = _node_num * 3 / 2;
2236
2237      if (!_matching) {
2238        _matching = new MatchingMap(_graph);
2239      }
2240      if (!_node_potential) {
2241        _node_potential = new NodePotential(_graph);
2242      }
2243      if (!_blossom_set) {
2244        _blossom_index = new IntNodeMap(_graph);
2245        _blossom_set = new BlossomSet(*_blossom_index);
2246        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2247      }
2248
2249      if (!_node_index) {
2250        _node_index = new IntNodeMap(_graph);
2251        _node_heap_index = new IntArcMap(_graph);
2252        _node_data = new RangeMap<NodeData>(_node_num,
2253                                            NodeData(*_node_heap_index));
2254      }
2255
2256      if (!_tree_set) {
2257        _tree_set_index = new IntIntMap(_blossom_num);
2258        _tree_set = new TreeSet(*_tree_set_index);
2259      }
2260      if (!_delta2) {
2261        _delta2_index = new IntIntMap(_blossom_num);
2262        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2263      }
2264      if (!_delta3) {
2265        _delta3_index = new IntEdgeMap(_graph);
2266        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
2267      }
2268      if (!_delta4) {
2269        _delta4_index = new IntIntMap(_blossom_num);
2270        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2271      }
2272    }
2273
2274    void destroyStructures() {
2275      if (_matching) {
2276        delete _matching;
2277      }
2278      if (_node_potential) {
2279        delete _node_potential;
2280      }
2281      if (_blossom_set) {
2282        delete _blossom_index;
2283        delete _blossom_set;
2284        delete _blossom_data;
2285      }
2286
2287      if (_node_index) {
2288        delete _node_index;
2289        delete _node_heap_index;
2290        delete _node_data;
2291      }
2292
2293      if (_tree_set) {
2294        delete _tree_set_index;
2295        delete _tree_set;
2296      }
2297      if (_delta2) {
2298        delete _delta2_index;
2299        delete _delta2;
2300      }
2301      if (_delta3) {
2302        delete _delta3_index;
2303        delete _delta3;
2304      }
2305      if (_delta4) {
2306        delete _delta4_index;
2307        delete _delta4;
2308      }
2309    }
2310
2311    void matchedToEven(int blossom, int tree) {
2312      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2313        _delta2->erase(blossom);
2314      }
2315
2316      if (!_blossom_set->trivial(blossom)) {
2317        (*_blossom_data)[blossom].pot -=
2318          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2319      }
2320
2321      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2322           n != INVALID; ++n) {
2323
2324        _blossom_set->increase(n, std::numeric_limits<Value>::max());
2325        int ni = (*_node_index)[n];
2326
2327        (*_node_data)[ni].heap.clear();
2328        (*_node_data)[ni].heap_index.clear();
2329
2330        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
2331
2332        for (InArcIt e(_graph, n); e != INVALID; ++e) {
2333          Node v = _graph.source(e);
2334          int vb = _blossom_set->find(v);
2335          int vi = (*_node_index)[v];
2336
2337          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2338            dualScale * _weight[e];
2339
2340          if ((*_blossom_data)[vb].status == EVEN) {
2341            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2342              _delta3->push(e, rw / 2);
2343            }
2344          } else {
2345            typename std::map<int, Arc>::iterator it =
2346              (*_node_data)[vi].heap_index.find(tree);
2347
2348            if (it != (*_node_data)[vi].heap_index.end()) {
2349              if ((*_node_data)[vi].heap[it->second] > rw) {
2350                (*_node_data)[vi].heap.replace(it->second, e);
2351                (*_node_data)[vi].heap.decrease(e, rw);
2352                it->second = e;
2353              }
2354            } else {
2355              (*_node_data)[vi].heap.push(e, rw);
2356              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2357            }
2358
2359            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2360              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2361
2362              if ((*_blossom_data)[vb].status == MATCHED) {
2363                if (_delta2->state(vb) != _delta2->IN_HEAP) {
2364                  _delta2->push(vb, _blossom_set->classPrio(vb) -
2365                               (*_blossom_data)[vb].offset);
2366                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2367                           (*_blossom_data)[vb].offset){
2368                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2369                                   (*_blossom_data)[vb].offset);
2370                }
2371              }
2372            }
2373          }
2374        }
2375      }
2376      (*_blossom_data)[blossom].offset = 0;
2377    }
2378
2379    void matchedToOdd(int blossom) {
2380      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2381        _delta2->erase(blossom);
2382      }
2383      (*_blossom_data)[blossom].offset += _delta_sum;
2384      if (!_blossom_set->trivial(blossom)) {
2385        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
2386                     (*_blossom_data)[blossom].offset);
2387      }
2388    }
2389
2390    void evenToMatched(int blossom, int tree) {
2391      if (!_blossom_set->trivial(blossom)) {
2392        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
2393      }
2394
2395      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2396           n != INVALID; ++n) {
2397        int ni = (*_node_index)[n];
2398        (*_node_data)[ni].pot -= _delta_sum;
2399
2400        for (InArcIt e(_graph, n); e != INVALID; ++e) {
2401          Node v = _graph.source(e);
2402          int vb = _blossom_set->find(v);
2403          int vi = (*_node_index)[v];
2404
2405          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2406            dualScale * _weight[e];
2407
2408          if (vb == blossom) {
2409            if (_delta3->state(e) == _delta3->IN_HEAP) {
2410              _delta3->erase(e);
2411            }
2412          } else if ((*_blossom_data)[vb].status == EVEN) {
2413
2414            if (_delta3->state(e) == _delta3->IN_HEAP) {
2415              _delta3->erase(e);
2416            }
2417
2418            int vt = _tree_set->find(vb);
2419
2420            if (vt != tree) {
2421
2422              Arc r = _graph.oppositeArc(e);
2423
2424              typename std::map<int, Arc>::iterator it =
2425                (*_node_data)[ni].heap_index.find(vt);
2426
2427              if (it != (*_node_data)[ni].heap_index.end()) {
2428                if ((*_node_data)[ni].heap[it->second] > rw) {
2429                  (*_node_data)[ni].heap.replace(it->second, r);
2430                  (*_node_data)[ni].heap.decrease(r, rw);
2431                  it->second = r;
2432                }
2433              } else {
2434                (*_node_data)[ni].heap.push(r, rw);
2435                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2436              }
2437
2438              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2439                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2440
2441                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2442                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2443                               (*_blossom_data)[blossom].offset);
2444                } else if ((*_delta2)[blossom] >
2445                           _blossom_set->classPrio(blossom) -
2446                           (*_blossom_data)[blossom].offset){
2447                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2448                                   (*_blossom_data)[blossom].offset);
2449                }
2450              }
2451            }
2452          } else {
2453
2454            typename std::map<int, Arc>::iterator it =
2455              (*_node_data)[vi].heap_index.find(tree);
2456
2457            if (it != (*_node_data)[vi].heap_index.end()) {
2458              (*_node_data)[vi].heap.erase(it->second);
2459              (*_node_data)[vi].heap_index.erase(it);
2460              if ((*_node_data)[vi].heap.empty()) {
2461                _blossom_set->increase(v, std::numeric_limits<Value>::max());
2462              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
2463                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
2464              }
2465
2466              if ((*_blossom_data)[vb].status == MATCHED) {
2467                if (_blossom_set->classPrio(vb) ==
2468                    std::numeric_limits<Value>::max()) {
2469                  _delta2->erase(vb);
2470                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
2471                           (*_blossom_data)[vb].offset) {
2472                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
2473                                   (*_blossom_data)[vb].offset);
2474                }
2475              }
2476            }
2477          }
2478        }
2479      }
2480    }
2481
2482    void oddToMatched(int blossom) {
2483      (*_blossom_data)[blossom].offset -= _delta_sum;
2484
2485      if (_blossom_set->classPrio(blossom) !=
2486          std::numeric_limits<Value>::max()) {
2487        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2488                       (*_blossom_data)[blossom].offset);
2489      }
2490
2491      if (!_blossom_set->trivial(blossom)) {
2492        _delta4->erase(blossom);
2493      }
2494    }
2495
2496    void oddToEven(int blossom, int tree) {
2497      if (!_blossom_set->trivial(blossom)) {
2498        _delta4->erase(blossom);
2499        (*_blossom_data)[blossom].pot -=
2500          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
2501      }
2502
2503      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2504           n != INVALID; ++n) {
2505        int ni = (*_node_index)[n];
2506
2507        _blossom_set->increase(n, std::numeric_limits<Value>::max());
2508
2509        (*_node_data)[ni].heap.clear();
2510        (*_node_data)[ni].heap_index.clear();
2511        (*_node_data)[ni].pot +=
2512          2 * _delta_sum - (*_blossom_data)[blossom].offset;
2513
2514        for (InArcIt e(_graph, n); e != INVALID; ++e) {
2515          Node v = _graph.source(e);
2516          int vb = _blossom_set->find(v);
2517          int vi = (*_node_index)[v];
2518
2519          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2520            dualScale * _weight[e];
2521
2522          if ((*_blossom_data)[vb].status == EVEN) {
2523            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2524              _delta3->push(e, rw / 2);
2525            }
2526          } else {
2527
2528            typename std::map<int, Arc>::iterator it =
2529              (*_node_data)[vi].heap_index.find(tree);
2530
2531            if (it != (*_node_data)[vi].heap_index.end()) {
2532              if ((*_node_data)[vi].heap[it->second] > rw) {
2533                (*_node_data)[vi].heap.replace(it->second, e);
2534                (*_node_data)[vi].heap.decrease(e, rw);
2535                it->second = e;
2536              }
2537            } else {
2538              (*_node_data)[vi].heap.push(e, rw);
2539              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2540            }
2541
2542            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2543              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2544
2545              if ((*_blossom_data)[vb].status == MATCHED) {
2546                if (_delta2->state(vb) != _delta2->IN_HEAP) {
2547                  _delta2->push(vb, _blossom_set->classPrio(vb) -
2548                               (*_blossom_data)[vb].offset);
2549                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2550                           (*_blossom_data)[vb].offset) {
2551                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2552                                   (*_blossom_data)[vb].offset);
2553                }
2554              }
2555            }
2556          }
2557        }
2558      }
2559      (*_blossom_data)[blossom].offset = 0;
2560    }
2561
2562    void alternatePath(int even, int tree) {
2563      int odd;
2564
2565      evenToMatched(even, tree);
2566      (*_blossom_data)[even].status = MATCHED;
2567
2568      while ((*_blossom_data)[even].pred != INVALID) {
2569        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
2570        (*_blossom_data)[odd].status = MATCHED;
2571        oddToMatched(odd);
2572        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2573
2574        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
2575        (*_blossom_data)[even].status = MATCHED;
2576        evenToMatched(even, tree);
2577        (*_blossom_data)[even].next =
2578          _graph.oppositeArc((*_blossom_data)[odd].pred);
2579      }
2580
2581    }
2582
2583    void destroyTree(int tree) {
2584      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2585        if ((*_blossom_data)[b].status == EVEN) {
2586          (*_blossom_data)[b].status = MATCHED;
2587          evenToMatched(b, tree);
2588        } else if ((*_blossom_data)[b].status == ODD) {
2589          (*_blossom_data)[b].status = MATCHED;
2590          oddToMatched(b);
2591        }
2592      }
2593      _tree_set->eraseClass(tree);
2594    }
2595
2596    void augmentOnEdge(const Edge& edge) {
2597
2598      int left = _blossom_set->find(_graph.u(edge));
2599      int right = _blossom_set->find(_graph.v(edge));
2600
2601      int left_tree = _tree_set->find(left);
2602      alternatePath(left, left_tree);
2603      destroyTree(left_tree);
2604
2605      int right_tree = _tree_set->find(right);
2606      alternatePath(right, right_tree);
2607      destroyTree(right_tree);
2608
2609      (*_blossom_data)[left].next = _graph.direct(edge, true);
2610      (*_blossom_data)[right].next = _graph.direct(edge, false);
2611    }
2612
2613    void extendOnArc(const Arc& arc) {
2614      int base = _blossom_set->find(_graph.target(arc));
2615      int tree = _tree_set->find(base);
2616
2617      int odd = _blossom_set->find(_graph.source(arc));
2618      _tree_set->insert(odd, tree);
2619      (*_blossom_data)[odd].status = ODD;
2620      matchedToOdd(odd);
2621      (*_blossom_data)[odd].pred = arc;
2622
2623      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
2624      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2625      _tree_set->insert(even, tree);
2626      (*_blossom_data)[even].status = EVEN;
2627      matchedToEven(even, tree);
2628    }
2629
2630    void shrinkOnEdge(const Edge& edge, int tree) {
2631      int nca = -1;
2632      std::vector<int> left_path, right_path;
2633
2634      {
2635        std::set<int> left_set, right_set;
2636        int left = _blossom_set->find(_graph.u(edge));
2637        left_path.push_back(left);
2638        left_set.insert(left);
2639
2640        int right = _blossom_set->find(_graph.v(edge));
2641        right_path.push_back(right);
2642        right_set.insert(right);
2643
2644        while (true) {
2645
2646          if ((*_blossom_data)[left].pred == INVALID) break;
2647
2648          left =
2649            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2650          left_path.push_back(left);
2651          left =
2652            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2653          left_path.push_back(left);
2654
2655          left_set.insert(left);
2656
2657          if (right_set.find(left) != right_set.end()) {
2658            nca = left;
2659            break;
2660          }
2661
2662          if ((*_blossom_data)[right].pred == INVALID) break;
2663
2664          right =
2665            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2666          right_path.push_back(right);
2667          right =
2668            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2669          right_path.push_back(right);
2670
2671          right_set.insert(right);
2672
2673          if (left_set.find(right) != left_set.end()) {
2674            nca = right;
2675            break;
2676          }
2677
2678        }
2679
2680        if (nca == -1) {
2681          if ((*_blossom_data)[left].pred == INVALID) {
2682            nca = right;
2683            while (left_set.find(nca) == left_set.end()) {
2684              nca =
2685                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2686              right_path.push_back(nca);
2687              nca =
2688                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2689              right_path.push_back(nca);
2690            }
2691          } else {
2692            nca = left;
2693            while (right_set.find(nca) == right_set.end()) {
2694              nca =
2695                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2696              left_path.push_back(nca);
2697              nca =
2698                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2699              left_path.push_back(nca);
2700            }
2701          }
2702        }
2703      }
2704
2705      std::vector<int> subblossoms;
2706      Arc prev;
2707
2708      prev = _graph.direct(edge, true);
2709      for (int i = 0; left_path[i] != nca; i += 2) {
2710        subblossoms.push_back(left_path[i]);
2711        (*_blossom_data)[left_path[i]].next = prev;
2712        _tree_set->erase(left_path[i]);
2713
2714        subblossoms.push_back(left_path[i + 1]);
2715        (*_blossom_data)[left_path[i + 1]].status = EVEN;
2716        oddToEven(left_path[i + 1], tree);
2717        _tree_set->erase(left_path[i + 1]);
2718        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
2719      }
2720
2721      int k = 0;
2722      while (right_path[k] != nca) ++k;
2723
2724      subblossoms.push_back(nca);
2725      (*_blossom_data)[nca].next = prev;
2726
2727      for (int i = k - 2; i >= 0; i -= 2) {
2728        subblossoms.push_back(right_path[i + 1]);
2729        (*_blossom_data)[right_path[i + 1]].status = EVEN;
2730        oddToEven(right_path[i + 1], tree);
2731        _tree_set->erase(right_path[i + 1]);
2732
2733        (*_blossom_data)[right_path[i + 1]].next =
2734          (*_blossom_data)[right_path[i + 1]].pred;
2735
2736        subblossoms.push_back(right_path[i]);
2737        _tree_set->erase(right_path[i]);
2738      }
2739
2740      int surface =
2741        _blossom_set->join(subblossoms.begin(), subblossoms.end());
2742
2743      for (int i = 0; i < int(subblossoms.size()); ++i) {
2744        if (!_blossom_set->trivial(subblossoms[i])) {
2745          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
2746        }
2747        (*_blossom_data)[subblossoms[i]].status = MATCHED;
2748      }
2749
2750      (*_blossom_data)[surface].pot = -2 * _delta_sum;
2751      (*_blossom_data)[surface].offset = 0;
2752      (*_blossom_data)[surface].status = EVEN;
2753      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
2754      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
2755
2756      _tree_set->insert(surface, tree);
2757      _tree_set->erase(nca);
2758    }
2759
2760    void splitBlossom(int blossom) {
2761      Arc next = (*_blossom_data)[blossom].next;
2762      Arc pred = (*_blossom_data)[blossom].pred;
2763
2764      int tree = _tree_set->find(blossom);
2765
2766      (*_blossom_data)[blossom].status = MATCHED;
2767      oddToMatched(blossom);
2768      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2769        _delta2->erase(blossom);
2770      }
2771
2772      std::vector<int> subblossoms;
2773      _blossom_set->split(blossom, std::back_inserter(subblossoms));
2774
2775      Value offset = (*_blossom_data)[blossom].offset;
2776      int b = _blossom_set->find(_graph.source(pred));
2777      int d = _blossom_set->find(_graph.source(next));
2778
2779      int ib = -1, id = -1;
2780      for (int i = 0; i < int(subblossoms.size()); ++i) {
2781        if (subblossoms[i] == b) ib = i;
2782        if (subblossoms[i] == d) id = i;
2783
2784        (*_blossom_data)[subblossoms[i]].offset = offset;
2785        if (!_blossom_set->trivial(subblossoms[i])) {
2786          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2787        }
2788        if (_blossom_set->classPrio(subblossoms[i]) !=
2789            std::numeric_limits<Value>::max()) {
2790          _delta2->push(subblossoms[i],
2791                        _blossom_set->classPrio(subblossoms[i]) -
2792                        (*_blossom_data)[subblossoms[i]].offset);
2793        }
2794      }
2795
2796      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2797        for (int i = (id + 1) % subblossoms.size();
2798             i != ib; i = (i + 2) % subblossoms.size()) {
2799          int sb = subblossoms[i];
2800          int tb = subblossoms[(i + 1) % subblossoms.size()];
2801          (*_blossom_data)[sb].next =
2802            _graph.oppositeArc((*_blossom_data)[tb].next);
2803        }
2804
2805        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2806          int sb = subblossoms[i];
2807          int tb = subblossoms[(i + 1) % subblossoms.size()];
2808          int ub = subblossoms[(i + 2) % subblossoms.size()];
2809
2810          (*_blossom_data)[sb].status = ODD;
2811          matchedToOdd(sb);
2812          _tree_set->insert(sb, tree);
2813          (*_blossom_data)[sb].pred = pred;
2814          (*_blossom_data)[sb].next =
2815                           _graph.oppositeArc((*_blossom_data)[tb].next);
2816
2817          pred = (*_blossom_data)[ub].next;
2818
2819          (*_blossom_data)[tb].status = EVEN;
2820          matchedToEven(tb, tree);
2821          _tree_set->insert(tb, tree);
2822          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2823        }
2824
2825        (*_blossom_data)[subblossoms[id]].status = ODD;
2826        matchedToOdd(subblossoms[id]);
2827        _tree_set->insert(subblossoms[id], tree);
2828        (*_blossom_data)[subblossoms[id]].next = next;
2829        (*_blossom_data)[subblossoms[id]].pred = pred;
2830
2831      } else {
2832
2833        for (int i = (ib + 1) % subblossoms.size();
2834             i != id; i = (i + 2) % subblossoms.size()) {
2835          int sb = subblossoms[i];
2836          int tb = subblossoms[(i + 1) % subblossoms.size()];
2837          (*_blossom_data)[sb].next =
2838            _graph.oppositeArc((*_blossom_data)[tb].next);
2839        }
2840
2841        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2842          int sb = subblossoms[i];
2843          int tb = subblossoms[(i + 1) % subblossoms.size()];
2844          int ub = subblossoms[(i + 2) % subblossoms.size()];
2845
2846          (*_blossom_data)[sb].status = ODD;
2847          matchedToOdd(sb);
2848          _tree_set->insert(sb, tree);
2849          (*_blossom_data)[sb].next = next;
2850          (*_blossom_data)[sb].pred =
2851            _graph.oppositeArc((*_blossom_data)[tb].next);
2852
2853          (*_blossom_data)[tb].status = EVEN;
2854          matchedToEven(tb, tree);
2855          _tree_set->insert(tb, tree);
2856          (*_blossom_data)[tb].pred =
2857            (*_blossom_data)[tb].next =
2858            _graph.oppositeArc((*_blossom_data)[ub].next);
2859          next = (*_blossom_data)[ub].next;
2860        }
2861
2862        (*_blossom_data)[subblossoms[ib]].status = ODD;
2863        matchedToOdd(subblossoms[ib]);
2864        _tree_set->insert(subblossoms[ib], tree);
2865        (*_blossom_data)[subblossoms[ib]].next = next;
2866        (*_blossom_data)[subblossoms[ib]].pred = pred;
2867      }
2868      _tree_set->erase(blossom);
2869    }
2870
2871    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2872      if (_blossom_set->trivial(blossom)) {
2873        int bi = (*_node_index)[base];
2874        Value pot = (*_node_data)[bi].pot;
2875
2876        (*_matching)[base] = matching;
2877        _blossom_node_list.push_back(base);
2878        (*_node_potential)[base] = pot;
2879      } else {
2880
2881        Value pot = (*_blossom_data)[blossom].pot;
2882        int bn = _blossom_node_list.size();
2883
2884        std::vector<int> subblossoms;
2885        _blossom_set->split(blossom, std::back_inserter(subblossoms));
2886        int b = _blossom_set->find(base);
2887        int ib = -1;
2888        for (int i = 0; i < int(subblossoms.size()); ++i) {
2889          if (subblossoms[i] == b) { ib = i; break; }
2890        }
2891
2892        for (int i = 1; i < int(subblossoms.size()); i += 2) {
2893          int sb = subblossoms[(ib + i) % subblossoms.size()];
2894          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2895
2896          Arc m = (*_blossom_data)[tb].next;
2897          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2898          extractBlossom(tb, _graph.source(m), m);
2899        }
2900        extractBlossom(subblossoms[ib], base, matching);
2901
2902        int en = _blossom_node_list.size();
2903
2904        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2905      }
2906    }
2907
2908    void extractMatching() {
2909      std::vector<int> blossoms;
2910      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2911        blossoms.push_back(c);
2912      }
2913
2914      for (int i = 0; i < int(blossoms.size()); ++i) {
2915
2916        Value offset = (*_blossom_data)[blossoms[i]].offset;
2917        (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2918        for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2919             n != INVALID; ++n) {
2920          (*_node_data)[(*_node_index)[n]].pot -= offset;
2921        }
2922
2923        Arc matching = (*_blossom_data)[blossoms[i]].next;
2924        Node base = _graph.source(matching);
2925        extractBlossom(blossoms[i], base, matching);
2926      }
2927    }
2928
2929  public:
2930
2931    /// \brief Constructor
2932    ///
2933    /// Constructor.
2934    MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2935      : _graph(graph), _weight(weight), _matching(0),
2936        _node_potential(0), _blossom_potential(), _blossom_node_list(),
2937        _node_num(0), _blossom_num(0),
2938
2939        _blossom_index(0), _blossom_set(0), _blossom_data(0),
2940        _node_index(0), _node_heap_index(0), _node_data(0),
2941        _tree_set_index(0), _tree_set(0),
2942
2943        _delta2_index(0), _delta2(0),
2944        _delta3_index(0), _delta3(0),
2945        _delta4_index(0), _delta4(0),
2946
2947        _delta_sum(), _unmatched(0),
2948
2949        _fractional(0)
2950    {}
2951
2952    ~MaxWeightedPerfectMatching() {
2953      destroyStructures();
2954      if (_fractional) {
2955        delete _fractional;
2956      }
2957    }
2958
2959    /// \name Execution Control
2960    /// The simplest way to execute the algorithm is to use the
2961    /// \ref run() member function.
2962
2963    ///@{
2964
2965    /// \brief Initialize the algorithm
2966    ///
2967    /// This function initializes the algorithm.
2968    void init() {
2969      createStructures();
2970
2971      for (ArcIt e(_graph); e != INVALID; ++e) {
2972        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
2973      }
2974      for (EdgeIt e(_graph); e != INVALID; ++e) {
2975        (*_delta3_index)[e] = _delta3->PRE_HEAP;
2976      }
2977      for (int i = 0; i < _blossom_num; ++i) {
2978        (*_delta2_index)[i] = _delta2->PRE_HEAP;
2979        (*_delta4_index)[i] = _delta4->PRE_HEAP;
2980      }
2981
2982      _unmatched = _node_num;
2983
2984      int index = 0;
2985      for (NodeIt n(_graph); n != INVALID; ++n) {
2986        Value max = - std::numeric_limits<Value>::max();
2987        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2988          if (_graph.target(e) == n) continue;
2989          if ((dualScale * _weight[e]) / 2 > max) {
2990            max = (dualScale * _weight[e]) / 2;
2991          }
2992        }
2993        (*_node_index)[n] = index;
2994        (*_node_data)[index].pot = max;
2995        int blossom =
2996          _blossom_set->insert(n, std::numeric_limits<Value>::max());
2997
2998        _tree_set->insert(blossom);
2999
3000        (*_blossom_data)[blossom].status = EVEN;
3001        (*_blossom_data)[blossom].pred = INVALID;
3002        (*_blossom_data)[blossom].next = INVALID;
3003        (*_blossom_data)[blossom].pot = 0;
3004        (*_blossom_data)[blossom].offset = 0;
3005        ++index;
3006      }
3007      for (EdgeIt e(_graph); e != INVALID; ++e) {
3008        int si = (*_node_index)[_graph.u(e)];
3009        int ti = (*_node_index)[_graph.v(e)];
3010        if (_graph.u(e) != _graph.v(e)) {
3011          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
3012                            dualScale * _weight[e]) / 2);
3013        }
3014      }
3015    }
3016
3017    /// \brief Initialize the algorithm with fractional matching
3018    ///
3019    /// This function initializes the algorithm with a fractional
3020    /// matching. This initialization is also called jumpstart heuristic.
3021    void fractionalInit() {
3022      createStructures();
3023     
3024      if (_fractional == 0) {
3025        _fractional = new FractionalMatching(_graph, _weight, false);
3026      }
3027      if (!_fractional->run()) {
3028        _unmatched = -1;
3029        return;
3030      }
3031
3032      for (ArcIt e(_graph); e != INVALID; ++e) {
3033        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
3034      }
3035      for (EdgeIt e(_graph); e != INVALID; ++e) {
3036        (*_delta3_index)[e] = _delta3->PRE_HEAP;
3037      }
3038      for (int i = 0; i < _blossom_num; ++i) {
3039        (*_delta2_index)[i] = _delta2->PRE_HEAP;
3040        (*_delta4_index)[i] = _delta4->PRE_HEAP;
3041      }
3042
3043      _unmatched = 0;
3044
3045      int index = 0;
3046      for (NodeIt n(_graph); n != INVALID; ++n) {
3047        Value pot = _fractional->nodeValue(n);
3048        (*_node_index)[n] = index;
3049        (*_node_data)[index].pot = pot;
3050        int blossom =
3051          _blossom_set->insert(n, std::numeric_limits<Value>::max());
3052
3053        (*_blossom_data)[blossom].status = MATCHED;
3054        (*_blossom_data)[blossom].pred = INVALID;
3055        (*_blossom_data)[blossom].next = _fractional->matching(n);
3056        (*_blossom_data)[blossom].pot = 0;
3057        (*_blossom_data)[blossom].offset = 0;
3058        ++index;
3059      }
3060
3061      typename Graph::template NodeMap<bool> processed(_graph, false);
3062      for (NodeIt n(_graph); n != INVALID; ++n) {
3063        if (processed[n]) continue;
3064        processed[n] = true;
3065        if (_fractional->matching(n) == INVALID) continue;
3066        int num = 1;
3067        Node v = _graph.target(_fractional->matching(n));
3068        while (n != v) {
3069          processed[v] = true;
3070          v = _graph.target(_fractional->matching(v));
3071          ++num;
3072        }
3073
3074        if (num % 2 == 1) {
3075          std::vector<int> subblossoms(num);
3076
3077          subblossoms[--num] = _blossom_set->find(n);
3078          v = _graph.target(_fractional->matching(n));
3079          while (n != v) {
3080            subblossoms[--num] = _blossom_set->find(v);
3081            v = _graph.target(_fractional->matching(v));           
3082          }
3083         
3084          int surface =
3085            _blossom_set->join(subblossoms.begin(), subblossoms.end());
3086          (*_blossom_data)[surface].status = EVEN;
3087          (*_blossom_data)[surface].pred = INVALID;
3088          (*_blossom_data)[surface].next = INVALID;
3089          (*_blossom_data)[surface].pot = 0;
3090          (*_blossom_data)[surface].offset = 0;
3091         
3092          _tree_set->insert(surface);
3093          ++_unmatched;
3094        }
3095      }
3096
3097      for (EdgeIt e(_graph); e != INVALID; ++e) {
3098        int si = (*_node_index)[_graph.u(e)];
3099        int sb = _blossom_set->find(_graph.u(e));
3100        int ti = (*_node_index)[_graph.v(e)];
3101        int tb = _blossom_set->find(_graph.v(e));
3102        if ((*_blossom_data)[sb].status == EVEN &&
3103            (*_blossom_data)[tb].status == EVEN && sb != tb) {
3104          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
3105                            dualScale * _weight[e]) / 2);
3106        }
3107      }
3108
3109      for (NodeIt n(_graph); n != INVALID; ++n) {
3110        int nb = _blossom_set->find(n);
3111        if ((*_blossom_data)[nb].status != MATCHED) continue;
3112        int ni = (*_node_index)[n];
3113
3114        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
3115          Node v = _graph.target(e);
3116          int vb = _blossom_set->find(v);
3117          int vi = (*_node_index)[v];
3118
3119          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
3120            dualScale * _weight[e];
3121
3122          if ((*_blossom_data)[vb].status == EVEN) {
3123
3124            int vt = _tree_set->find(vb);
3125
3126            typename std::map<int, Arc>::iterator it =
3127              (*_node_data)[ni].heap_index.find(vt);
3128
3129            if (it != (*_node_data)[ni].heap_index.end()) {
3130              if ((*_node_data)[ni].heap[it->second] > rw) {
3131                (*_node_data)[ni].heap.replace(it->second, e);
3132                (*_node_data)[ni].heap.decrease(e, rw);
3133                it->second = e;
3134              }
3135            } else {
3136              (*_node_data)[ni].heap.push(e, rw);
3137              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
3138            }
3139          }
3140        }
3141           
3142        if (!(*_node_data)[ni].heap.empty()) {
3143          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
3144          _delta2->push(nb, _blossom_set->classPrio(nb));
3145        }
3146      }
3147    }
3148
3149    /// \brief Start the algorithm
3150    ///
3151    /// This function starts the algorithm.
3152    ///
3153    /// \pre \ref init() or \ref fractionalInit() must be called before
3154    /// using this function.
3155    bool start() {
3156      enum OpType {
3157        D2, D3, D4
3158      };
3159
3160      if (_unmatched == -1) return false;
3161
3162      while (_unmatched > 0) {
3163        Value d2 = !_delta2->empty() ?
3164          _delta2->prio() : std::numeric_limits<Value>::max();
3165
3166        Value d3 = !_delta3->empty() ?
3167          _delta3->prio() : std::numeric_limits<Value>::max();
3168
3169        Value d4 = !_delta4->empty() ?
3170          _delta4->prio() : std::numeric_limits<Value>::max();
3171
3172        _delta_sum = d3; OpType ot = D3;
3173        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
3174        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
3175
3176        if (_delta_sum == std::numeric_limits<Value>::max()) {
3177          return false;
3178        }
3179
3180        switch (ot) {
3181        case D2:
3182          {
3183            int blossom = _delta2->top();
3184            Node n = _blossom_set->classTop(blossom);
3185            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
3186            extendOnArc(e);
3187          }
3188          break;
3189        case D3:
3190          {
3191            Edge e = _delta3->top();
3192
3193            int left_blossom = _blossom_set->find(_graph.u(e));
3194            int right_blossom = _blossom_set->find(_graph.v(e));
3195
3196            if (left_blossom == right_blossom) {
3197              _delta3->pop();
3198            } else {
3199              int left_tree = _tree_set->find(left_blossom);
3200              int right_tree = _tree_set->find(right_blossom);
3201
3202              if (left_tree == right_tree) {
3203                shrinkOnEdge(e, left_tree);
3204              } else {
3205                augmentOnEdge(e);
3206                _unmatched -= 2;
3207              }
3208            }
3209          } break;
3210        case D4:
3211          splitBlossom(_delta4->top());
3212          break;
3213        }
3214      }
3215      extractMatching();
3216      return true;
3217    }
3218
3219    /// \brief Run the algorithm.
3220    ///
3221    /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
3222    ///
3223    /// \note mwpm.run() is just a shortcut of the following code.
3224    /// \code
3225    ///   mwpm.fractionalInit();
3226    ///   mwpm.start();
3227    /// \endcode
3228    bool run() {
3229      fractionalInit();
3230      return start();
3231    }
3232
3233    /// @}
3234
3235    /// \name Primal Solution
3236    /// Functions to get the primal solution, i.e. the maximum weighted
3237    /// perfect matching.\n
3238    /// Either \ref run() or \ref start() function should be called before
3239    /// using them.
3240
3241    /// @{
3242
3243    /// \brief Return the weight of the matching.
3244    ///
3245    /// This function returns the weight of the found matching.
3246    ///
3247    /// \pre Either run() or start() must be called before using this function.
3248    Value matchingWeight() const {
3249      Value sum = 0;
3250      for (NodeIt n(_graph); n != INVALID; ++n) {
3251        if ((*_matching)[n] != INVALID) {
3252          sum += _weight[(*_matching)[n]];
3253        }
3254      }
3255      return sum / 2;
3256    }
3257
3258    /// \brief Return \c true if the given edge is in the matching.
3259    ///
3260    /// This function returns \c true if the given edge is in the found
3261    /// matching.
3262    ///
3263    /// \pre Either run() or start() must be called before using this function.
3264    bool matching(const Edge& edge) const {
3265      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
3266    }
3267
3268    /// \brief Return the matching arc (or edge) incident to the given node.
3269    ///
3270    /// This function returns the matching arc (or edge) incident to the
3271    /// given node in the found matching or \c INVALID if the node is
3272    /// not covered by the matching.
3273    ///
3274    /// \pre Either run() or start() must be called before using this function.
3275    Arc matching(const Node& node) const {
3276      return (*_matching)[node];
3277    }
3278
3279    /// \brief Return a const reference to the matching map.
3280    ///
3281    /// This function returns a const reference to a node map that stores
3282    /// the matching arc (or edge) incident to each node.
3283    const MatchingMap& matchingMap() const {
3284      return *_matching;
3285    }
3286
3287    /// \brief Return the mate of the given node.
3288    ///
3289    /// This function returns the mate of the given node in the found
3290    /// matching or \c INVALID if the node is not covered by the matching.
3291    ///
3292    /// \pre Either run() or start() must be called before using this function.
3293    Node mate(const Node& node) const {
3294      return _graph.target((*_matching)[node]);
3295    }
3296
3297    /// @}
3298
3299    /// \name Dual Solution
3300    /// Functions to get the dual solution.\n
3301    /// Either \ref run() or \ref start() function should be called before
3302    /// using them.
3303
3304    /// @{
3305
3306    /// \brief Return the value of the dual solution.
3307    ///
3308    /// This function returns the value of the dual solution.
3309    /// It should be equal to the primal value scaled by \ref dualScale
3310    /// "dual scale".
3311    ///
3312    /// \pre Either run() or start() must be called before using this function.
3313    Value dualValue() const {
3314      Value sum = 0;
3315      for (NodeIt n(_graph); n != INVALID; ++n) {
3316        sum += nodeValue(n);
3317      }
3318      for (int i = 0; i < blossomNum(); ++i) {
3319        sum += blossomValue(i) * (blossomSize(i) / 2);
3320      }
3321      return sum;
3322    }
3323
3324    /// \brief Return the dual value (potential) of the given node.
3325    ///
3326    /// This function returns the dual value (potential) of the given node.
3327    ///
3328    /// \pre Either run() or start() must be called before using this function.
3329    Value nodeValue(const Node& n) const {
3330      return (*_node_potential)[n];
3331    }
3332
3333    /// \brief Return the number of the blossoms in the basis.
3334    ///
3335    /// This function returns the number of the blossoms in the basis.
3336    ///
3337    /// \pre Either run() or start() must be called before using this function.
3338    /// \see BlossomIt
3339    int blossomNum() const {
3340      return _blossom_potential.size();
3341    }
3342
3343    /// \brief Return the number of the nodes in the given blossom.
3344    ///
3345    /// This function returns the number of the nodes in the given blossom.
3346    ///
3347    /// \pre Either run() or start() must be called before using this function.
3348    /// \see BlossomIt
3349    int blossomSize(int k) const {
3350      return _blossom_potential[k].end - _blossom_potential[k].begin;
3351    }
3352
3353    /// \brief Return the dual value (ptential) of the given blossom.
3354    ///
3355    /// This function returns the dual value (ptential) of the given blossom.
3356    ///
3357    /// \pre Either run() or start() must be called before using this function.
3358    Value blossomValue(int k) const {
3359      return _blossom_potential[k].value;
3360    }
3361
3362    /// \brief Iterator for obtaining the nodes of a blossom.
3363    ///
3364    /// This class provides an iterator for obtaining the nodes of the
3365    /// given blossom. It lists a subset of the nodes.
3366    /// Before using this iterator, you must allocate a
3367    /// MaxWeightedPerfectMatching class and execute it.
3368    class BlossomIt {
3369    public:
3370
3371      /// \brief Constructor.
3372      ///
3373      /// Constructor to get the nodes of the given variable.
3374      ///
3375      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
3376      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
3377      /// must be called before initializing this iterator.
3378      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
3379        : _algorithm(&algorithm)
3380      {
3381        _index = _algorithm->_blossom_potential[variable].begin;
3382        _last = _algorithm->_blossom_potential[variable].end;
3383      }
3384
3385      /// \brief Conversion to \c Node.
3386      ///
3387      /// Conversion to \c Node.
3388      operator Node() const {
3389        return _algorithm->_blossom_node_list[_index];
3390      }
3391
3392      /// \brief Increment operator.
3393      ///
3394      /// Increment operator.
3395      BlossomIt& operator++() {
3396        ++_index;
3397        return *this;
3398      }
3399
3400      /// \brief Validity checking
3401      ///
3402      /// This function checks whether the iterator is invalid.
3403      bool operator==(Invalid) const { return _index == _last; }
3404
3405      /// \brief Validity checking
3406      ///
3407      /// This function checks whether the iterator is valid.
3408      bool operator!=(Invalid) const { return _index != _last; }
3409
3410    private:
3411      const MaxWeightedPerfectMatching* _algorithm;
3412      int _last;
3413      int _index;
3414    };
3415
3416    /// @}
3417
3418  };
3419
3420} //END OF NAMESPACE LEMON
3421
3422#endif //LEMON_MATCHING_H
Note: See TracBrowser for help on using the repository browser.