COIN-OR::LEMON - Graph Library

Ticket #48: 64ad48007fb2.patch

File 64ad48007fb2.patch, 114.6 KB (added by Balazs Dezso, 16 years ago)

Port of the matching algorithms

  • lemon/Makefile.am

    # HG changeset patch
    # User Balazs Dezso <deba@inf.elte.hu>
    # Date 1223898960 -7200
    # Node ID 64ad48007fb28697c78f3f4c644894cd9572f14d
    # Parent  6dbd5184c6a99e936a7de899e64c114e30c0f94f
    Port maximum matching algorithms from svn 3498 (ticket #48)
    
    diff -r 6dbd5184c6a9 -r 64ad48007fb2 lemon/Makefile.am
    a b  
    3535        lemon/list_graph.h \
    3636        lemon/maps.h \
    3737        lemon/math.h \
     38        lemon/max_matching.h \
    3839        lemon/path.h \
    3940        lemon/random.h \
    4041        lemon/smart_graph.h \
  • new file lemon/max_matching.h

    diff -r 6dbd5184c6a9 -r 64ad48007fb2 lemon/max_matching.h
    - +  
     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
  • test/CMakeLists.txt

    diff -r 6dbd5184c6a9 -r 64ad48007fb2 test/CMakeLists.txt
    a b  
    1616  heap_test
    1717  kruskal_test
    1818  maps_test
     19  max_matching_test
     20  max_weighted_matching_test
    1921  random_test
    2022  path_test
    2123  time_measure_test
  • test/Makefile.am

    diff -r 6dbd5184c6a9 -r 64ad48007fb2 test/Makefile.am
    a b  
    1919        test/heap_test \
    2020        test/kruskal_test \
    2121        test/maps_test \
     22        test/max_matching_test \
     23        test/max_weighted_matching_test \
    2224        test/random_test \
    2325        test/path_test \
    2426        test/test_tools_fail \
     
    4244test_heap_test_SOURCES = test/heap_test.cc
    4345test_kruskal_test_SOURCES = test/kruskal_test.cc
    4446test_maps_test_SOURCES = test/maps_test.cc
     47test_max_matching_test_SOURCES = test/max_matching_test.cc
     48test_max_weighted_matching_test_SOURCES = test/max_weighted_matching_test.cc
    4549test_path_test_SOURCES = test/path_test.cc
    4650test_random_test_SOURCES = test/random_test.cc
    4751test_test_tools_fail_SOURCES = test/test_tools_fail.cc
  • new file test/max_matching_test.cc

    diff -r 6dbd5184c6a9 -r 64ad48007fb2 test/max_matching_test.cc
    - +  
     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#include <iostream>
     20#include <vector>
     21#include <queue>
     22#include <lemon/math.h>
     23#include <cstdlib>
     24
     25#include "test_tools.h"
     26#include <lemon/list_graph.h>
     27#include <lemon/max_matching.h>
     28
     29using namespace std;
     30using namespace lemon;
     31
     32int main() {
     33
     34  typedef ListGraph Graph;
     35
     36  GRAPH_TYPEDEFS(Graph);
     37
     38  Graph g;
     39  g.clear();
     40
     41  std::vector<Graph::Node> nodes;
     42  for (int i=0; i<13; ++i)
     43      nodes.push_back(g.addNode());
     44
     45  g.addEdge(nodes[0], nodes[0]);
     46  g.addEdge(nodes[6], nodes[10]);
     47  g.addEdge(nodes[5], nodes[10]);
     48  g.addEdge(nodes[4], nodes[10]);
     49  g.addEdge(nodes[3], nodes[11]);
     50  g.addEdge(nodes[1], nodes[6]);
     51  g.addEdge(nodes[4], nodes[7]);
     52  g.addEdge(nodes[1], nodes[8]);
     53  g.addEdge(nodes[0], nodes[8]);
     54  g.addEdge(nodes[3], nodes[12]);
     55  g.addEdge(nodes[6], nodes[9]);
     56  g.addEdge(nodes[9], nodes[11]);
     57  g.addEdge(nodes[2], nodes[10]);
     58  g.addEdge(nodes[10], nodes[8]);
     59  g.addEdge(nodes[5], nodes[8]);
     60  g.addEdge(nodes[6], nodes[3]);
     61  g.addEdge(nodes[0], nodes[5]);
     62  g.addEdge(nodes[6], nodes[12]);
     63
     64  MaxMatching<Graph> max_matching(g);
     65  max_matching.init();
     66  max_matching.startDense();
     67
     68  int s=0;
     69  Graph::NodeMap<Node> mate(g,INVALID);
     70  max_matching.mateMap(mate);
     71  for(NodeIt v(g); v!=INVALID; ++v) {
     72    if ( mate[v]!=INVALID ) ++s;
     73  }
     74  int size=int(s/2);  //size will be used as the size of a maxmatching
     75
     76  for(NodeIt v(g); v!=INVALID; ++v) {
     77    max_matching.mate(v);
     78  }
     79
     80  check ( size == max_matching.size(), "mate() returns a different size matching than max_matching.size()" );
     81
     82  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos0(g);
     83  max_matching.decomposition(pos0);
     84
     85  max_matching.init();
     86  max_matching.startSparse();
     87  s=0;
     88  max_matching.mateMap(mate);
     89  for(NodeIt v(g); v!=INVALID; ++v) {
     90    if ( mate[v]!=INVALID ) ++s;
     91  }
     92  check ( int(s/2) == size, "The size does not equal!" );
     93
     94  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos1(g);
     95  max_matching.decomposition(pos1);
     96
     97  max_matching.run();
     98  s=0;
     99  max_matching.mateMap(mate);
     100  for(NodeIt v(g); v!=INVALID; ++v) {
     101    if ( mate[v]!=INVALID ) ++s;
     102  }
     103  check ( int(s/2) == size, "The size does not equal!" );
     104
     105  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos2(g);
     106  max_matching.decomposition(pos2);
     107
     108  max_matching.run();
     109  s=0;
     110  max_matching.mateMap(mate);
     111  for(NodeIt v(g); v!=INVALID; ++v) {
     112    if ( mate[v]!=INVALID ) ++s;
     113  }
     114  check ( int(s/2) == size, "The size does not equal!" );
     115
     116  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos(g);
     117  max_matching.decomposition(pos);
     118
     119  bool ismatching=true;
     120  for(NodeIt v(g); v!=INVALID; ++v) {
     121    if ( mate[v]!=INVALID ) {
     122      Node u=mate[v];
     123      if (mate[u]!=v) ismatching=false;
     124    }
     125  }
     126  check ( ismatching, "It is not a matching!" );
     127
     128  bool coincide=true;
     129  for(NodeIt v(g); v!=INVALID; ++v) {
     130   if ( pos0[v] != pos1[v] || pos1[v]!=pos2[v] || pos2[v]!=pos[v] ) {
     131     coincide=false;
     132    }
     133  }
     134  check ( coincide, "The decompositions do not coincide! " );
     135
     136  bool noarc=true;
     137  for(EdgeIt e(g); e!=INVALID; ++e) {
     138   if ( (pos[g.v(e)]==max_matching.C &&
     139         pos[g.u(e)]==max_matching.D) ||
     140         (pos[g.v(e)]==max_matching.D &&
     141          pos[g.u(e)]==max_matching.C) )
     142      noarc=false;
     143  }
     144  check ( noarc, "There are arcs between D and C!" );
     145
     146  bool oddcomp=true;
     147  Graph::NodeMap<bool> todo(g,true);
     148  int num_comp=0;
     149  for(NodeIt v(g); v!=INVALID; ++v) {
     150   if ( pos[v]==max_matching.D && todo[v] ) {
     151      int comp_size=1;
     152      ++num_comp;
     153      std::queue<Node> Q;
     154      Q.push(v);
     155      todo.set(v,false);
     156      while (!Q.empty()) {
     157        Node w=Q.front();
     158        Q.pop();
     159        for(IncEdgeIt e(g,w); e!=INVALID; ++e) {
     160          Node u=g.runningNode(e);
     161          if ( pos[u]==max_matching.D && todo[u] ) {
     162            ++comp_size;
     163            Q.push(u);
     164            todo.set(u,false);
     165          }
     166        }
     167      }
     168      if ( !(comp_size % 2) ) oddcomp=false;
     169    }
     170  }
     171  check ( oddcomp, "A component of g[D] is not odd." );
     172
     173  int barrier=0;
     174  for(NodeIt v(g); v!=INVALID; ++v) {
     175    if ( pos[v]==max_matching.A ) ++barrier;
     176  }
     177  int expected_size=int( countNodes(g)-num_comp+barrier)/2;
     178  check ( size==expected_size, "The size of the matching is wrong." );
     179
     180  return 0;
     181}
  • new file test/max_weighted_matching_test.cc

    diff -r 6dbd5184c6a9 -r 64ad48007fb2 test/max_weighted_matching_test.cc
    - +  
     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#include <iostream>
     20#include <sstream>
     21#include <vector>
     22#include <queue>
     23#include <lemon/math.h>
     24#include <cstdlib>
     25
     26#include "test_tools.h"
     27#include <lemon/smart_graph.h>
     28#include <lemon/max_matching.h>
     29#include <lemon/lgf_reader.h>
     30
     31using namespace std;
     32using namespace lemon;
     33
     34GRAPH_TYPEDEFS(SmartGraph);
     35
     36const int lgfn = 3;
     37const std::string lgf[lgfn] = {
     38  "@nodes\n"
     39  "label\n"
     40  "0\n"
     41  "1\n"
     42  "2\n"
     43  "3\n"
     44  "4\n"
     45  "5\n"
     46  "6\n"
     47  "7\n"
     48  "@edges\n"
     49  "label        weight\n"
     50  "7    4       0       984\n"
     51  "0    7       1       73\n"
     52  "7    1       2       204\n"
     53  "2    3       3       583\n"
     54  "2    7       4       565\n"
     55  "2    1       5       582\n"
     56  "0    4       6       551\n"
     57  "2    5       7       385\n"
     58  "1    5       8       561\n"
     59  "5    3       9       484\n"
     60  "7    5       10      904\n"
     61  "3    6       11      47\n"
     62  "7    6       12      888\n"
     63  "3    0       13      747\n"
     64  "6    1       14      310\n",
     65
     66  "@nodes\n"
     67  "label\n"
     68  "0\n"
     69  "1\n"
     70  "2\n"
     71  "3\n"
     72  "4\n"
     73  "5\n"
     74  "6\n"
     75  "7\n"
     76  "@edges\n"
     77  "label        weight\n"
     78  "2    5       0       710\n"
     79  "0    5       1       241\n"
     80  "2    4       2       856\n"
     81  "2    6       3       762\n"
     82  "4    1       4       747\n"
     83  "6    1       5       962\n"
     84  "4    7       6       723\n"
     85  "1    7       7       661\n"
     86  "2    3       8       376\n"
     87  "1    0       9       416\n"
     88  "6    7       10      391\n",
     89
     90  "@nodes\n"
     91  "label\n"
     92  "0\n"
     93  "1\n"
     94  "2\n"
     95  "3\n"
     96  "4\n"
     97  "5\n"
     98  "6\n"
     99  "7\n"
     100  "@edges\n"
     101  "label        weight\n"
     102  "6    2       0       553\n"
     103  "0    7       1       653\n"
     104  "6    3       2       22\n"
     105  "4    7       3       846\n"
     106  "7    2       4       981\n"
     107  "7    6       5       250\n"
     108  "5    2       6       539\n"
     109};
     110
     111void checkMatching(const SmartGraph& graph,
     112                   const SmartGraph::EdgeMap<int>& weight,
     113                   const MaxWeightedMatching<SmartGraph>& mwm) {
     114  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
     115    if (graph.u(e) == graph.v(e)) continue;
     116    int rw =
     117      mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
     118
     119    for (int i = 0; i < mwm.blossomNum(); ++i) {
     120      bool s = false, t = false;
     121      for (MaxWeightedMatching<SmartGraph>::BlossomIt n(mwm, i);
     122           n != INVALID; ++n) {
     123        if (graph.u(e) == n) s = true;
     124        if (graph.v(e) == n) t = true;
     125      }
     126      if (s == true && t == true) {
     127        rw += mwm.blossomValue(i);
     128      }
     129    }
     130    rw -= weight[e] * mwm.dualScale;
     131
     132    check(rw >= 0, "Negative reduced weight");
     133    check(rw == 0 || !mwm.matching(e),
     134          "Non-zero reduced weight on matching arc");
     135  }
     136
     137  int pv = 0;
     138  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     139    if (mwm.matching(n) != INVALID) {
     140      check(mwm.nodeValue(n) >= 0, "Invalid node value");
     141      pv += weight[mwm.matching(n)];
     142      SmartGraph::Node o = graph.target(mwm.matching(n));
     143      check(mwm.mate(n) == o, "Invalid matching");
     144      check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
     145            "Invalid matching");
     146    } else {
     147      check(mwm.mate(n) == INVALID, "Invalid matching");
     148      check(mwm.nodeValue(n) == 0, "Invalid matching");
     149    }
     150  }
     151
     152  int dv = 0;
     153  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     154    dv += mwm.nodeValue(n);
     155  }
     156
     157  for (int i = 0; i < mwm.blossomNum(); ++i) {
     158    check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
     159    check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
     160    dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
     161  }
     162
     163  check(pv * mwm.dualScale == dv * 2, "Wrong duality");
     164
     165  return;
     166}
     167
     168void checkPerfectMatching(const SmartGraph& graph,
     169                          const SmartGraph::EdgeMap<int>& weight,
     170                          const MaxWeightedPerfectMatching<SmartGraph>& mwpm) {
     171  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
     172    if (graph.u(e) == graph.v(e)) continue;
     173    int rw =
     174      mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
     175
     176    for (int i = 0; i < mwpm.blossomNum(); ++i) {
     177      bool s = false, t = false;
     178      for (MaxWeightedPerfectMatching<SmartGraph>::BlossomIt n(mwpm, i);
     179           n != INVALID; ++n) {
     180        if (graph.u(e) == n) s = true;
     181        if (graph.v(e) == n) t = true;
     182      }
     183      if (s == true && t == true) {
     184        rw += mwpm.blossomValue(i);
     185      }
     186    }
     187    rw -= weight[e] * mwpm.dualScale;
     188
     189    check(rw >= 0, "Negative reduced weight");
     190    check(rw == 0 || !mwpm.matching(e),
     191          "Non-zero reduced weight on matching arc");
     192  }
     193
     194  int pv = 0;
     195  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     196    check(mwpm.matching(n) != INVALID, "Non perfect");
     197    pv += weight[mwpm.matching(n)];
     198    SmartGraph::Node o = graph.target(mwpm.matching(n));
     199    check(mwpm.mate(n) == o, "Invalid matching");
     200    check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
     201          "Invalid matching");
     202  }
     203
     204  int dv = 0;
     205  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     206    dv += mwpm.nodeValue(n);
     207  }
     208
     209  for (int i = 0; i < mwpm.blossomNum(); ++i) {
     210    check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
     211    check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
     212    dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
     213  }
     214
     215  check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
     216
     217  return;
     218}
     219
     220
     221int main() {
     222
     223  for (int i = 0; i < lgfn; ++i) {
     224    SmartGraph graph;
     225    SmartGraph::EdgeMap<int> weight(graph);
     226
     227    istringstream lgfs(lgf[i]);
     228    graphReader(graph, lgfs).
     229      edgeMap("weight", weight).run();
     230
     231    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
     232    mwm.run();
     233    checkMatching(graph, weight, mwm);
     234
     235    MaxMatching<SmartGraph> mm(graph);
     236    mm.run();
     237
     238    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
     239    bool perfect = mwpm.run();
     240
     241    check(perfect == (mm.size() * 2 == countNodes(graph)),
     242          "Perfect matching found");
     243
     244    if (perfect) {
     245      checkPerfectMatching(graph, weight, mwpm);
     246    }
     247  }
     248
     249  return 0;
     250}