COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/max_matching.h @ 2553:bfced05fa852

Last change on this file since 2553:bfced05fa852 was 2553:bfced05fa852, checked in by Alpar Juttner, 11 years ago

Happy New Year to LEMON (+ better update-copyright-header script)

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