COIN-OR::LEMON - Graph Library

source: lemon-main/lemon/max_matching.h @ 326:64ad48007fb2

Last change on this file since 326:64ad48007fb2 was 326:64ad48007fb2, checked in by Balazs Dezso <deba@…>, 16 years ago

Port maximum matching algorithms from svn 3498 (ticket #48)

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