COIN-OR::LEMON - Graph Library

Changeset 2548:a3ba22ebccc6 in lemon-0.x


Ignore:
Timestamp:
12/28/07 12:00:51 (11 years ago)
Author:
Balazs Dezso
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3425
Message:

Edmond's Blossom shrinking algroithm:
MaxWeightedMatching?
MaxWeightedPerfectMatching?

Files:
4 edited

Legend:

Unmodified
Added
Removed
  • doc/groups.dox

    r2530 r2548  
    319319for find matchings in graphs and bipartite graphs.
    320320
    321 This group provides some algorithm objects and function
    322 to calculate matchings in graphs and bipartite graphs.
     321This group provides some algorithm objects and function to calculate
     322matchings in graphs and bipartite graphs. The general matching problem is
     323finding a subset of the edges which does not shares common endpoints.
     324 
     325There are several different algorithms for calculate matchings in
     326graphs.  The matching problems in bipartite graphs are generally
     327easier than in general graphs. The goal of the matching optimization
     328can be the finding maximum cardinality, maximum weight or minimum cost
     329matching. The search can be constrained to find perfect or
     330maximum cardinality matching.
     331
     332Lemon contains the next algorithms:
     333- \ref lemon::MaxBipartiteMatching "MaxBipartiteMatching" Hopcroft-Karp
     334  augmenting path algorithm for calculate maximum cardinality matching in
     335  bipartite graphs
     336- \ref lemon::PrBipartiteMatching "PrBipartiteMatching" Push-Relabel
     337  algorithm for calculate maximum cardinality matching in bipartite graphs
     338- \ref lemon::MaxWeightedBipartiteMatching "MaxWeightedBipartiteMatching"
     339  Successive shortest path algorithm for calculate maximum weighted matching
     340  and maximum weighted bipartite matching in bipartite graph
     341- \ref lemon::MinCostMaxBipartiteMatching "MinCostMaxBipartiteMatching"
     342  Successive shortest path algorithm for calculate minimum cost maximum
     343  matching in bipartite graph
     344- \ref lemon::MaxMatching "MaxMatching" Edmond's blossom shrinking algorithm
     345  for calculate maximum cardinality matching in general graph
     346- \ref lemon::MaxWeightedMatching "MaxWeightedMatching" Edmond's blossom
     347  shrinking algorithm for calculate maximum weighted matching in general
     348  graph
     349- \ref lemon::MaxWeightedPerfectMatching "MaxWeightedPerfectMatching"
     350  Edmond's blossom shrinking algorithm for calculate maximum weighted
     351  perfect matching in general graph
    323352
    324353\image html bipartite_matching.png
  • lemon/bin_heap.h

    r2547 r2548  
    5353
    5454  public:
     55    ///\e
    5556    typedef _ItemIntMap ItemIntMap;
     57    ///\e
    5658    typedef _Prio Prio;
     59    ///\e
    5760    typedef typename ItemIntMap::Key Item;
     61    ///\e
    5862    typedef std::pair<Item,Prio> Pair;
     63    ///\e
    5964    typedef _Compare Compare;
    6065
     
    322327    }
    323328
     329    /// \brief Replaces an item in the heap.
     330    ///
     331    /// The \c i item is replaced with \c j item. The \c i item should
     332    /// be in the heap, while the \c j should be out of the heap. The
     333    /// \c i item will out of the heap and \c j will be in the heap
     334    /// with the same prioriority as prevoiusly the \c i item.
     335    void replace(const Item& i, const Item& j) {
     336      int idx = iim[i];
     337      iim.set(i, iim[j]);
     338      iim.set(j, idx);
     339      data[idx].first = j;
     340    }
     341
    324342  }; // class BinHeap
    325343 
  • lemon/max_matching.h

    r2505 r2548  
    2020#define LEMON_MAX_MATCHING_H
    2121
     22#include <vector>
    2223#include <queue>
     24#include <set>
     25
    2326#include <lemon/bits/invalid.h>
    2427#include <lemon/unionfind.h>
    2528#include <lemon/graph_utils.h>
     29#include <lemon/bin_heap.h>
    2630
    2731///\ingroup matching
     
    3236
    3337  ///\ingroup matching
    34 
     38  ///
    3539  ///\brief Edmonds' alternating forest maximum matching algorithm.
    3640  ///
     
    619623
    620624  };
     625
     626  /// \ingroup matching
     627  ///
     628  /// \brief Weighted matching in general undirected graphs
     629  ///
     630  /// This class provides an efficient implementation of Edmond's
     631  /// maximum weighted matching algorithm. The implementation is based
     632  /// on extensive use of priority queues and provides
     633  /// \f$O(nm\log(n))\f$ time complexity.
     634  ///
     635  /// The maximum weighted matching problem is to find undirected
     636  /// edges in the graph with maximum overall weight and no two of
     637  /// them shares their endpoints. The problem can be formulated with
     638  /// the next linear program:
     639  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
     640  ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
     641  /// \f[x_e \ge 0\quad \forall e\in E\f]
     642  /// \f[\max \sum_{e\in E}x_ew_e\f]
     643  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
     644  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in
     645  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
     646  /// the nodes.
     647  ///
     648  /// The algorithm calculates an optimal matching and a proof of the
     649  /// optimality. The solution of the dual problem can be used to check
     650  /// the result of the algorithm. The dual linear problem is the next:
     651  /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
     652  /// \f[y_u \ge 0 \quad \forall u \in V\f]
     653  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
     654  /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
     655  ///
     656  /// The algorithm can be executed with \c run() or the \c init() and
     657  /// then the \c start() member functions. After it the matching can
     658  /// be asked with \c matching() or mate() functions. The dual
     659  /// solution can be get with \c nodeValue(), \c blossomNum() and \c
     660  /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
     661  /// "BlossomIt" nested class which is able to iterate on the nodes
     662  /// of a blossom. If the value type is integral then the dual
     663  /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
     664  ///
     665  /// \author Balazs Dezso
     666  template <typename _UGraph,
     667            typename _WeightMap = typename _UGraph::template UEdgeMap<int> >
     668  class MaxWeightedMatching {
     669  public:
     670
     671    typedef _UGraph UGraph;
     672    typedef _WeightMap WeightMap;
     673    typedef typename WeightMap::Value Value;
     674
     675    /// \brief Scaling factor for dual solution
     676    ///
     677    /// Scaling factor for dual solution, it is equal to 4 or 1
     678    /// according to the value type.
     679    static const int dualScale =
     680      std::numeric_limits<Value>::is_integer ? 4 : 1;
     681
     682    typedef typename UGraph::template NodeMap<typename UGraph::Edge>
     683    MatchingMap;
     684   
     685  private:
     686
     687    UGRAPH_TYPEDEFS(typename UGraph);
     688
     689    typedef typename UGraph::template NodeMap<Value> NodePotential;
     690    typedef std::vector<Node> BlossomNodeList;
     691
     692    struct BlossomVariable {
     693      int begin, end;
     694      Value value;
     695     
     696      BlossomVariable(int _begin, int _end, Value _value)
     697        : begin(_begin), end(_end), value(_value) {}
     698
     699    };
     700
     701    typedef std::vector<BlossomVariable> BlossomPotential;
     702
     703    const UGraph& _ugraph;
     704    const WeightMap& _weight;
     705
     706    MatchingMap* _matching;
     707
     708    NodePotential* _node_potential;
     709
     710    BlossomPotential _blossom_potential;
     711    BlossomNodeList _blossom_node_list;
     712
     713    int _node_num;
     714    int _blossom_num;
     715
     716    typedef typename UGraph::template NodeMap<int> NodeIntMap;
     717    typedef typename UGraph::template EdgeMap<int> EdgeIntMap;
     718    typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap;
     719    typedef IntegerMap<int> IntIntMap;
     720
     721    enum Status {
     722      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
     723    };
     724
     725    typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
     726    struct BlossomData {
     727      int tree;
     728      Status status;
     729      Edge pred, next;
     730      Value pot, offset;
     731      Node base;
     732    };
     733
     734    NodeIntMap *_blossom_index;
     735    BlossomSet *_blossom_set;
     736    IntegerMap<BlossomData>* _blossom_data;
     737
     738    NodeIntMap *_node_index;
     739    EdgeIntMap *_node_heap_index;
     740
     741    struct NodeData {
     742     
     743      NodeData(EdgeIntMap& node_heap_index)
     744        : heap(node_heap_index) {}
     745     
     746      int blossom;
     747      Value pot;
     748      BinHeap<Value, EdgeIntMap> heap;
     749      std::map<int, Edge> heap_index;
     750     
     751      int tree;
     752    };
     753
     754    IntegerMap<NodeData>* _node_data;
     755
     756    typedef ExtendFindEnum<IntIntMap> TreeSet;
     757
     758    IntIntMap *_tree_set_index;
     759    TreeSet *_tree_set;
     760
     761    NodeIntMap *_delta1_index;
     762    BinHeap<Value, NodeIntMap> *_delta1;
     763
     764    IntIntMap *_delta2_index;
     765    BinHeap<Value, IntIntMap> *_delta2;
     766   
     767    UEdgeIntMap *_delta3_index;
     768    BinHeap<Value, UEdgeIntMap> *_delta3;
     769
     770    IntIntMap *_delta4_index;
     771    BinHeap<Value, IntIntMap> *_delta4;
     772
     773    Value _delta_sum;
     774
     775    void createStructures() {
     776      _node_num = countNodes(_ugraph);
     777      _blossom_num = _node_num * 3 / 2;
     778
     779      if (!_matching) {
     780        _matching = new MatchingMap(_ugraph);
     781      }
     782      if (!_node_potential) {
     783        _node_potential = new NodePotential(_ugraph);
     784      }
     785      if (!_blossom_set) {
     786        _blossom_index = new NodeIntMap(_ugraph);
     787        _blossom_set = new BlossomSet(*_blossom_index);
     788        _blossom_data = new IntegerMap<BlossomData>(_blossom_num);
     789      }
     790
     791      if (!_node_index) {
     792        _node_index = new NodeIntMap(_ugraph);
     793        _node_heap_index = new EdgeIntMap(_ugraph);
     794        _node_data = new IntegerMap<NodeData>(_node_num,
     795                                              NodeData(*_node_heap_index));
     796      }
     797
     798      if (!_tree_set) {
     799        _tree_set_index = new IntIntMap(_blossom_num);
     800        _tree_set = new TreeSet(*_tree_set_index);
     801      }
     802      if (!_delta1) {
     803        _delta1_index = new NodeIntMap(_ugraph);
     804        _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);
     805      }
     806      if (!_delta2) {
     807        _delta2_index = new IntIntMap(_blossom_num);
     808        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
     809      }
     810      if (!_delta3) {
     811        _delta3_index = new UEdgeIntMap(_ugraph);
     812        _delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index);
     813      }
     814      if (!_delta4) {
     815        _delta4_index = new IntIntMap(_blossom_num);
     816        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
     817      }
     818    }
     819
     820    void destroyStructures() {
     821      _node_num = countNodes(_ugraph);
     822      _blossom_num = _node_num * 3 / 2;
     823
     824      if (_matching) {
     825        delete _matching;
     826      }
     827      if (_node_potential) {
     828        delete _node_potential;
     829      }
     830      if (_blossom_set) {
     831        delete _blossom_index;
     832        delete _blossom_set;
     833        delete _blossom_data;
     834      }
     835
     836      if (_node_index) {
     837        delete _node_index;
     838        delete _node_heap_index;
     839        delete _node_data;                           
     840      }
     841
     842      if (_tree_set) {
     843        delete _tree_set_index;
     844        delete _tree_set;
     845      }
     846      if (_delta1) {
     847        delete _delta1_index;
     848        delete _delta1;
     849      }
     850      if (_delta2) {
     851        delete _delta2_index;
     852        delete _delta2;
     853      }
     854      if (_delta3) {
     855        delete _delta3_index;
     856        delete _delta3;
     857      }
     858      if (_delta4) {
     859        delete _delta4_index;
     860        delete _delta4;
     861      }
     862    }
     863
     864    void matchedToEven(int blossom, int tree) {
     865      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
     866        _delta2->erase(blossom);
     867      }
     868
     869      if (!_blossom_set->trivial(blossom)) {
     870        (*_blossom_data)[blossom].pot -=
     871          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
     872      }
     873
     874      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
     875           n != INVALID; ++n) {
     876
     877        _blossom_set->increase(n, std::numeric_limits<Value>::max());
     878        int ni = (*_node_index)[n];
     879
     880        (*_node_data)[ni].heap.clear();
     881        (*_node_data)[ni].heap_index.clear();
     882
     883        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
     884
     885        _delta1->push(n, (*_node_data)[ni].pot);
     886
     887        for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
     888          Node v = _ugraph.source(e);
     889          int vb = _blossom_set->find(v);
     890          int vi = (*_node_index)[v];
     891
     892          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
     893            dualScale * _weight[e];
     894
     895          if ((*_blossom_data)[vb].status == EVEN) {
     896            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
     897              _delta3->push(e, rw / 2);
     898            }
     899          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
     900            if (_delta3->state(e) != _delta3->IN_HEAP) {
     901              _delta3->push(e, rw);
     902            }       
     903          } else {
     904            typename std::map<int, Edge>::iterator it =
     905              (*_node_data)[vi].heap_index.find(tree);   
     906
     907            if (it != (*_node_data)[vi].heap_index.end()) {
     908              if ((*_node_data)[vi].heap[it->second] > rw) {
     909                (*_node_data)[vi].heap.replace(it->second, e);
     910                (*_node_data)[vi].heap.decrease(e, rw);
     911                it->second = e;
     912              }
     913            } else {
     914              (*_node_data)[vi].heap.push(e, rw);
     915              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
     916            }
     917
     918            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
     919              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
     920
     921              if ((*_blossom_data)[vb].status == MATCHED) {
     922                if (_delta2->state(vb) != _delta2->IN_HEAP) {
     923                  _delta2->push(vb, _blossom_set->classPrio(vb) -
     924                               (*_blossom_data)[vb].offset);
     925                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
     926                           (*_blossom_data)[vb].offset){
     927                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
     928                                   (*_blossom_data)[vb].offset);
     929                }
     930              }
     931            }
     932          }
     933        }
     934      }
     935      (*_blossom_data)[blossom].offset = 0;
     936    }
     937
     938    void matchedToOdd(int blossom) {
     939      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
     940        _delta2->erase(blossom);
     941      }
     942      (*_blossom_data)[blossom].offset += _delta_sum;
     943      if (!_blossom_set->trivial(blossom)) {
     944        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
     945                     (*_blossom_data)[blossom].offset);
     946      }
     947    }
     948
     949    void evenToMatched(int blossom, int tree) {
     950      if (!_blossom_set->trivial(blossom)) {
     951        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
     952      }
     953
     954      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
     955           n != INVALID; ++n) {
     956        int ni = (*_node_index)[n];
     957        (*_node_data)[ni].pot -= _delta_sum;
     958
     959        _delta1->erase(n);
     960
     961        for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
     962          Node v = _ugraph.source(e);
     963          int vb = _blossom_set->find(v);
     964          int vi = (*_node_index)[v];
     965
     966          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
     967            dualScale * _weight[e];
     968
     969          if (vb == blossom) {
     970            if (_delta3->state(e) == _delta3->IN_HEAP) {
     971              _delta3->erase(e);
     972            }
     973          } else if ((*_blossom_data)[vb].status == EVEN) {
     974
     975            if (_delta3->state(e) == _delta3->IN_HEAP) {
     976              _delta3->erase(e);
     977            }
     978
     979            int vt = _tree_set->find(vb);
     980           
     981            if (vt != tree) {
     982
     983              Edge r = _ugraph.oppositeEdge(e);
     984
     985              typename std::map<int, Edge>::iterator it =
     986                (*_node_data)[ni].heap_index.find(vt);   
     987
     988              if (it != (*_node_data)[ni].heap_index.end()) {
     989                if ((*_node_data)[ni].heap[it->second] > rw) {
     990                  (*_node_data)[ni].heap.replace(it->second, r);
     991                  (*_node_data)[ni].heap.decrease(r, rw);
     992                  it->second = r;
     993                }
     994              } else {
     995                (*_node_data)[ni].heap.push(r, rw);
     996                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
     997              }
     998             
     999              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
     1000                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
     1001               
     1002                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
     1003                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
     1004                               (*_blossom_data)[blossom].offset);
     1005                } else if ((*_delta2)[blossom] >
     1006                           _blossom_set->classPrio(blossom) -
     1007                           (*_blossom_data)[blossom].offset){
     1008                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
     1009                                   (*_blossom_data)[blossom].offset);
     1010                }
     1011              }
     1012            }
     1013           
     1014          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
     1015            if (_delta3->state(e) == _delta3->IN_HEAP) {
     1016              _delta3->erase(e);
     1017            }
     1018          } else {
     1019         
     1020            typename std::map<int, Edge>::iterator it =
     1021              (*_node_data)[vi].heap_index.find(tree);
     1022
     1023            if (it != (*_node_data)[vi].heap_index.end()) {
     1024              (*_node_data)[vi].heap.erase(it->second);
     1025              (*_node_data)[vi].heap_index.erase(it);
     1026              if ((*_node_data)[vi].heap.empty()) {
     1027                _blossom_set->increase(v, std::numeric_limits<Value>::max());
     1028              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
     1029                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
     1030              }
     1031             
     1032              if ((*_blossom_data)[vb].status == MATCHED) {
     1033                if (_blossom_set->classPrio(vb) ==
     1034                    std::numeric_limits<Value>::max()) {
     1035                  _delta2->erase(vb);
     1036                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
     1037                           (*_blossom_data)[vb].offset) {
     1038                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
     1039                                   (*_blossom_data)[vb].offset);
     1040                }
     1041              }
     1042            }   
     1043          }
     1044        }
     1045      }
     1046    }
     1047
     1048    void oddToMatched(int blossom) {
     1049      (*_blossom_data)[blossom].offset -= _delta_sum;
     1050
     1051      if (_blossom_set->classPrio(blossom) !=
     1052          std::numeric_limits<Value>::max()) {
     1053        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
     1054                       (*_blossom_data)[blossom].offset);
     1055      }
     1056
     1057      if (!_blossom_set->trivial(blossom)) {
     1058        _delta4->erase(blossom);
     1059      }
     1060    }
     1061
     1062    void oddToEven(int blossom, int tree) {
     1063      if (!_blossom_set->trivial(blossom)) {
     1064        _delta4->erase(blossom);
     1065        (*_blossom_data)[blossom].pot -=
     1066          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
     1067      }
     1068
     1069      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
     1070           n != INVALID; ++n) {
     1071        int ni = (*_node_index)[n];
     1072
     1073        _blossom_set->increase(n, std::numeric_limits<Value>::max());
     1074
     1075        (*_node_data)[ni].heap.clear();
     1076        (*_node_data)[ni].heap_index.clear();
     1077        (*_node_data)[ni].pot +=
     1078          2 * _delta_sum - (*_blossom_data)[blossom].offset;
     1079
     1080        _delta1->push(n, (*_node_data)[ni].pot);
     1081
     1082        for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
     1083          Node v = _ugraph.source(e);
     1084          int vb = _blossom_set->find(v);
     1085          int vi = (*_node_index)[v];
     1086
     1087          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
     1088            dualScale * _weight[e];
     1089
     1090          if ((*_blossom_data)[vb].status == EVEN) {
     1091            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
     1092              _delta3->push(e, rw / 2);
     1093            }
     1094          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
     1095            if (_delta3->state(e) != _delta3->IN_HEAP) {
     1096              _delta3->push(e, rw);
     1097            }
     1098          } else {
     1099         
     1100            typename std::map<int, Edge>::iterator it =
     1101              (*_node_data)[vi].heap_index.find(tree);   
     1102
     1103            if (it != (*_node_data)[vi].heap_index.end()) {
     1104              if ((*_node_data)[vi].heap[it->second] > rw) {
     1105                (*_node_data)[vi].heap.replace(it->second, e);
     1106                (*_node_data)[vi].heap.decrease(e, rw);
     1107                it->second = e;
     1108              }
     1109            } else {
     1110              (*_node_data)[vi].heap.push(e, rw);
     1111              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
     1112            }
     1113
     1114            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
     1115              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
     1116
     1117              if ((*_blossom_data)[vb].status == MATCHED) {
     1118                if (_delta2->state(vb) != _delta2->IN_HEAP) {
     1119                  _delta2->push(vb, _blossom_set->classPrio(vb) -
     1120                               (*_blossom_data)[vb].offset);
     1121                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
     1122                           (*_blossom_data)[vb].offset) {
     1123                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
     1124                                   (*_blossom_data)[vb].offset);
     1125                }
     1126              }
     1127            }
     1128          }
     1129        }
     1130      }
     1131      (*_blossom_data)[blossom].offset = 0;
     1132    }
     1133
     1134
     1135    void matchedToUnmatched(int blossom) {
     1136      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
     1137        _delta2->erase(blossom);
     1138      }
     1139
     1140      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
     1141           n != INVALID; ++n) {
     1142        int ni = (*_node_index)[n];
     1143
     1144        _blossom_set->increase(n, std::numeric_limits<Value>::max());
     1145
     1146        (*_node_data)[ni].heap.clear();
     1147        (*_node_data)[ni].heap_index.clear();
     1148
     1149        for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
     1150          Node v = _ugraph.target(e);
     1151          int vb = _blossom_set->find(v);
     1152          int vi = (*_node_index)[v];
     1153
     1154          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
     1155            dualScale * _weight[e];
     1156
     1157          if ((*_blossom_data)[vb].status == EVEN) {
     1158            if (_delta3->state(e) != _delta3->IN_HEAP) {
     1159              _delta3->push(e, rw);
     1160            }
     1161          }
     1162        }
     1163      }
     1164    }
     1165
     1166    void unmatchedToMatched(int blossom) {
     1167      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
     1168           n != INVALID; ++n) {
     1169        int ni = (*_node_index)[n];
     1170
     1171        for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
     1172          Node v = _ugraph.source(e);
     1173          int vb = _blossom_set->find(v);
     1174          int vi = (*_node_index)[v];
     1175
     1176          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
     1177            dualScale * _weight[e];
     1178
     1179          if (vb == blossom) {
     1180            if (_delta3->state(e) == _delta3->IN_HEAP) {
     1181              _delta3->erase(e);
     1182            }
     1183          } else if ((*_blossom_data)[vb].status == EVEN) {
     1184
     1185            if (_delta3->state(e) == _delta3->IN_HEAP) {
     1186              _delta3->erase(e);
     1187            }
     1188
     1189            int vt = _tree_set->find(vb);
     1190           
     1191            Edge r = _ugraph.oppositeEdge(e);
     1192           
     1193            typename std::map<int, Edge>::iterator it =
     1194              (*_node_data)[ni].heap_index.find(vt);     
     1195           
     1196            if (it != (*_node_data)[ni].heap_index.end()) {
     1197              if ((*_node_data)[ni].heap[it->second] > rw) {
     1198                (*_node_data)[ni].heap.replace(it->second, r);
     1199                (*_node_data)[ni].heap.decrease(r, rw);
     1200                it->second = r;
     1201              }
     1202            } else {
     1203              (*_node_data)[ni].heap.push(r, rw);
     1204              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
     1205            }
     1206             
     1207            if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
     1208              _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
     1209             
     1210              if (_delta2->state(blossom) != _delta2->IN_HEAP) {
     1211                _delta2->push(blossom, _blossom_set->classPrio(blossom) -
     1212                             (*_blossom_data)[blossom].offset);
     1213              } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
     1214                         (*_blossom_data)[blossom].offset){
     1215                _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
     1216                                 (*_blossom_data)[blossom].offset);
     1217              }
     1218            }
     1219           
     1220          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
     1221            if (_delta3->state(e) == _delta3->IN_HEAP) {
     1222              _delta3->erase(e);
     1223            }
     1224          }
     1225        }
     1226      }
     1227    }
     1228
     1229    void alternatePath(int even, int tree) {
     1230      int odd;
     1231
     1232      evenToMatched(even, tree);
     1233      (*_blossom_data)[even].status = MATCHED;
     1234
     1235      while ((*_blossom_data)[even].pred != INVALID) {
     1236        odd = _blossom_set->find(_ugraph.target((*_blossom_data)[even].pred));
     1237        (*_blossom_data)[odd].status = MATCHED;
     1238        oddToMatched(odd);
     1239        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
     1240     
     1241        even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred));
     1242        (*_blossom_data)[even].status = MATCHED;
     1243        evenToMatched(even, tree);
     1244        (*_blossom_data)[even].next =
     1245          _ugraph.oppositeEdge((*_blossom_data)[odd].pred);
     1246      }
     1247     
     1248    }
     1249
     1250    void destroyTree(int tree) {
     1251      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
     1252        if ((*_blossom_data)[b].status == EVEN) {
     1253          (*_blossom_data)[b].status = MATCHED;
     1254          evenToMatched(b, tree);
     1255        } else if ((*_blossom_data)[b].status == ODD) {
     1256          (*_blossom_data)[b].status = MATCHED;
     1257          oddToMatched(b);
     1258        }
     1259      }
     1260      _tree_set->eraseClass(tree);
     1261    }
     1262   
     1263
     1264    void unmatchNode(const Node& node) {
     1265      int blossom = _blossom_set->find(node);
     1266      int tree = _tree_set->find(blossom);
     1267
     1268      alternatePath(blossom, tree);
     1269      destroyTree(tree);
     1270
     1271      (*_blossom_data)[blossom].status = UNMATCHED;
     1272      (*_blossom_data)[blossom].base = node;
     1273      matchedToUnmatched(blossom);
     1274    }
     1275
     1276
     1277    void augmentOnEdge(const UEdge& edge) {
     1278     
     1279      int left = _blossom_set->find(_ugraph.source(edge));
     1280      int right = _blossom_set->find(_ugraph.target(edge));
     1281
     1282      if ((*_blossom_data)[left].status == EVEN) {
     1283        int left_tree = _tree_set->find(left);
     1284        alternatePath(left, left_tree);
     1285        destroyTree(left_tree);
     1286      } else {
     1287        (*_blossom_data)[left].status = MATCHED;
     1288        unmatchedToMatched(left);
     1289      }
     1290
     1291      if ((*_blossom_data)[right].status == EVEN) {
     1292        int right_tree = _tree_set->find(right);
     1293        alternatePath(right, right_tree);
     1294        destroyTree(right_tree);
     1295      } else {
     1296        (*_blossom_data)[right].status = MATCHED;
     1297        unmatchedToMatched(right);
     1298      }
     1299
     1300      (*_blossom_data)[left].next = _ugraph.direct(edge, true);
     1301      (*_blossom_data)[right].next = _ugraph.direct(edge, false);
     1302    }
     1303
     1304    void extendOnEdge(const Edge& edge) {
     1305      int base = _blossom_set->find(_ugraph.target(edge));
     1306      int tree = _tree_set->find(base);
     1307     
     1308      int odd = _blossom_set->find(_ugraph.source(edge));
     1309      _tree_set->insert(odd, tree);
     1310      (*_blossom_data)[odd].status = ODD;
     1311      matchedToOdd(odd);
     1312      (*_blossom_data)[odd].pred = edge;
     1313
     1314      int even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].next));
     1315      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
     1316      _tree_set->insert(even, tree);
     1317      (*_blossom_data)[even].status = EVEN;
     1318      matchedToEven(even, tree);
     1319    }
     1320   
     1321    void shrinkOnEdge(const UEdge& uedge, int tree) {
     1322      int nca = -1;
     1323      std::vector<int> left_path, right_path;
     1324
     1325      {
     1326        std::set<int> left_set, right_set;
     1327        int left = _blossom_set->find(_ugraph.source(uedge));
     1328        left_path.push_back(left);
     1329        left_set.insert(left);
     1330
     1331        int right = _blossom_set->find(_ugraph.target(uedge));
     1332        right_path.push_back(right);
     1333        right_set.insert(right);
     1334
     1335        while (true) {
     1336
     1337          if ((*_blossom_data)[left].pred == INVALID) break;
     1338
     1339          left =
     1340            _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
     1341          left_path.push_back(left);
     1342          left =
     1343            _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
     1344          left_path.push_back(left);
     1345
     1346          left_set.insert(left);
     1347
     1348          if (right_set.find(left) != right_set.end()) {
     1349            nca = left;
     1350            break;
     1351          }
     1352
     1353          if ((*_blossom_data)[right].pred == INVALID) break;
     1354
     1355          right =
     1356            _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
     1357          right_path.push_back(right);
     1358          right =
     1359            _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
     1360          right_path.push_back(right);
     1361
     1362          right_set.insert(right);
     1363 
     1364          if (left_set.find(right) != left_set.end()) {
     1365            nca = right;
     1366            break;
     1367          }
     1368
     1369        }
     1370
     1371        if (nca == -1) {
     1372          if ((*_blossom_data)[left].pred == INVALID) {
     1373            nca = right;
     1374            while (left_set.find(nca) == left_set.end()) {
     1375              nca =
     1376                _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
     1377              right_path.push_back(nca);
     1378              nca =
     1379                _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
     1380              right_path.push_back(nca);
     1381            }
     1382          } else {
     1383            nca = left;
     1384            while (right_set.find(nca) == right_set.end()) {
     1385              nca =
     1386                _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
     1387              left_path.push_back(nca);
     1388              nca =
     1389                _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
     1390              left_path.push_back(nca);
     1391            }
     1392          }
     1393        }
     1394      }
     1395
     1396      std::vector<int> subblossoms;
     1397      Edge prev;
     1398
     1399      prev = _ugraph.direct(uedge, true);
     1400      for (int i = 0; left_path[i] != nca; i += 2) {
     1401        subblossoms.push_back(left_path[i]);
     1402        (*_blossom_data)[left_path[i]].next = prev;
     1403        _tree_set->erase(left_path[i]);
     1404
     1405        subblossoms.push_back(left_path[i + 1]);
     1406        (*_blossom_data)[left_path[i + 1]].status = EVEN;
     1407        oddToEven(left_path[i + 1], tree);
     1408        _tree_set->erase(left_path[i + 1]);
     1409        prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred);
     1410      }
     1411
     1412      int k = 0;
     1413      while (right_path[k] != nca) ++k;
     1414
     1415      subblossoms.push_back(nca);
     1416      (*_blossom_data)[nca].next = prev;
     1417     
     1418      for (int i = k - 2; i >= 0; i -= 2) {
     1419        subblossoms.push_back(right_path[i + 1]);
     1420        (*_blossom_data)[right_path[i + 1]].status = EVEN;
     1421        oddToEven(right_path[i + 1], tree);
     1422        _tree_set->erase(right_path[i + 1]);
     1423
     1424        (*_blossom_data)[right_path[i + 1]].next =
     1425          (*_blossom_data)[right_path[i + 1]].pred;
     1426
     1427        subblossoms.push_back(right_path[i]);
     1428        _tree_set->erase(right_path[i]);
     1429      }
     1430
     1431      int surface =
     1432        _blossom_set->join(subblossoms.begin(), subblossoms.end());
     1433
     1434      for (int i = 0; i < int(subblossoms.size()); ++i) {
     1435        if (!_blossom_set->trivial(subblossoms[i])) {
     1436          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
     1437        }
     1438        (*_blossom_data)[subblossoms[i]].status = MATCHED;
     1439      }
     1440
     1441      (*_blossom_data)[surface].pot = -2 * _delta_sum;
     1442      (*_blossom_data)[surface].offset = 0;
     1443      (*_blossom_data)[surface].status = EVEN;
     1444      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
     1445      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
     1446
     1447      _tree_set->insert(surface, tree);
     1448      _tree_set->erase(nca);
     1449    }
     1450
     1451    void splitBlossom(int blossom) {
     1452      Edge next = (*_blossom_data)[blossom].next;
     1453      Edge pred = (*_blossom_data)[blossom].pred;
     1454
     1455      int tree = _tree_set->find(blossom);
     1456
     1457      (*_blossom_data)[blossom].status = MATCHED;
     1458      oddToMatched(blossom);
     1459      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
     1460        _delta2->erase(blossom);
     1461      }
     1462
     1463      std::vector<int> subblossoms;
     1464      _blossom_set->split(blossom, std::back_inserter(subblossoms));
     1465
     1466      Value offset = (*_blossom_data)[blossom].offset;
     1467      int b = _blossom_set->find(_ugraph.source(pred));
     1468      int d = _blossom_set->find(_ugraph.source(next));
     1469     
     1470      int ib, id;
     1471      for (int i = 0; i < int(subblossoms.size()); ++i) {
     1472        if (subblossoms[i] == b) ib = i;
     1473        if (subblossoms[i] == d) id = i;
     1474
     1475        (*_blossom_data)[subblossoms[i]].offset = offset;
     1476        if (!_blossom_set->trivial(subblossoms[i])) {
     1477          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
     1478        }
     1479        if (_blossom_set->classPrio(subblossoms[i]) !=
     1480            std::numeric_limits<Value>::max()) {
     1481          _delta2->push(subblossoms[i],
     1482                        _blossom_set->classPrio(subblossoms[i]) -
     1483                        (*_blossom_data)[subblossoms[i]].offset);
     1484        }
     1485      }
     1486     
     1487      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
     1488        for (int i = (id + 1) % subblossoms.size();
     1489             i != ib; i = (i + 2) % subblossoms.size()) {
     1490          int sb = subblossoms[i];
     1491          int tb = subblossoms[(i + 1) % subblossoms.size()];
     1492          (*_blossom_data)[sb].next =
     1493            _ugraph.oppositeEdge((*_blossom_data)[tb].next);
     1494        }
     1495
     1496        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
     1497          int sb = subblossoms[i];
     1498          int tb = subblossoms[(i + 1) % subblossoms.size()];
     1499          int ub = subblossoms[(i + 2) % subblossoms.size()];
     1500
     1501          (*_blossom_data)[sb].status = ODD;
     1502          matchedToOdd(sb);
     1503          _tree_set->insert(sb, tree);
     1504          (*_blossom_data)[sb].pred = pred;
     1505          (*_blossom_data)[sb].next =
     1506                           _ugraph.oppositeEdge((*_blossom_data)[tb].next);
     1507         
     1508          pred = (*_blossom_data)[ub].next;
     1509     
     1510          (*_blossom_data)[tb].status = EVEN;
     1511          matchedToEven(tb, tree);
     1512          _tree_set->insert(tb, tree);
     1513          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
     1514        }
     1515     
     1516        (*_blossom_data)[subblossoms[id]].status = ODD;
     1517        matchedToOdd(subblossoms[id]);
     1518        _tree_set->insert(subblossoms[id], tree);
     1519        (*_blossom_data)[subblossoms[id]].next = next;
     1520        (*_blossom_data)[subblossoms[id]].pred = pred;
     1521     
     1522      } else {
     1523
     1524        for (int i = (ib + 1) % subblossoms.size();
     1525             i != id; i = (i + 2) % subblossoms.size()) {
     1526          int sb = subblossoms[i];
     1527          int tb = subblossoms[(i + 1) % subblossoms.size()];
     1528          (*_blossom_data)[sb].next =
     1529            _ugraph.oppositeEdge((*_blossom_data)[tb].next);
     1530        }       
     1531
     1532        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
     1533          int sb = subblossoms[i];
     1534          int tb = subblossoms[(i + 1) % subblossoms.size()];
     1535          int ub = subblossoms[(i + 2) % subblossoms.size()];
     1536
     1537          (*_blossom_data)[sb].status = ODD;
     1538          matchedToOdd(sb);
     1539          _tree_set->insert(sb, tree);
     1540          (*_blossom_data)[sb].next = next;
     1541          (*_blossom_data)[sb].pred =
     1542            _ugraph.oppositeEdge((*_blossom_data)[tb].next);
     1543
     1544          (*_blossom_data)[tb].status = EVEN;
     1545          matchedToEven(tb, tree);
     1546          _tree_set->insert(tb, tree);
     1547          (*_blossom_data)[tb].pred =
     1548            (*_blossom_data)[tb].next =
     1549            _ugraph.oppositeEdge((*_blossom_data)[ub].next);
     1550          next = (*_blossom_data)[ub].next;
     1551        }
     1552
     1553        (*_blossom_data)[subblossoms[ib]].status = ODD;
     1554        matchedToOdd(subblossoms[ib]);
     1555        _tree_set->insert(subblossoms[ib], tree);
     1556        (*_blossom_data)[subblossoms[ib]].next = next;
     1557        (*_blossom_data)[subblossoms[ib]].pred = pred;
     1558      }
     1559      _tree_set->erase(blossom);
     1560    }
     1561
     1562    void extractBlossom(int blossom, const Node& base, const Edge& matching) {
     1563      if (_blossom_set->trivial(blossom)) {
     1564        int bi = (*_node_index)[base];
     1565        Value pot = (*_node_data)[bi].pot;
     1566
     1567        _matching->set(base, matching);
     1568        _blossom_node_list.push_back(base);
     1569        _node_potential->set(base, pot);
     1570      } else {
     1571
     1572        Value pot = (*_blossom_data)[blossom].pot;
     1573        int bn = _blossom_node_list.size();
     1574       
     1575        std::vector<int> subblossoms;
     1576        _blossom_set->split(blossom, std::back_inserter(subblossoms));
     1577        int b = _blossom_set->find(base);
     1578        int ib = -1;
     1579        for (int i = 0; i < int(subblossoms.size()); ++i) {
     1580          if (subblossoms[i] == b) { ib = i; break; }
     1581        }
     1582       
     1583        for (int i = 1; i < int(subblossoms.size()); i += 2) {
     1584          int sb = subblossoms[(ib + i) % subblossoms.size()];
     1585          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
     1586         
     1587          Edge m = (*_blossom_data)[tb].next;
     1588          extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m));
     1589          extractBlossom(tb, _ugraph.source(m), m);
     1590        }
     1591        extractBlossom(subblossoms[ib], base, matching);     
     1592       
     1593        int en = _blossom_node_list.size();
     1594       
     1595        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
     1596      }
     1597    }
     1598
     1599    void extractMatching() {
     1600      std::vector<int> blossoms;
     1601      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
     1602        blossoms.push_back(c);
     1603      }
     1604
     1605      for (int i = 0; i < int(blossoms.size()); ++i) {
     1606        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
     1607
     1608          Value offset = (*_blossom_data)[blossoms[i]].offset;
     1609          (*_blossom_data)[blossoms[i]].pot += 2 * offset;
     1610          for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
     1611               n != INVALID; ++n) {
     1612            (*_node_data)[(*_node_index)[n]].pot -= offset;
     1613          }
     1614
     1615          Edge matching = (*_blossom_data)[blossoms[i]].next;
     1616          Node base = _ugraph.source(matching);
     1617          extractBlossom(blossoms[i], base, matching);
     1618        } else {
     1619          Node base = (*_blossom_data)[blossoms[i]].base;         
     1620          extractBlossom(blossoms[i], base, INVALID);
     1621        }
     1622      }
     1623    }
     1624   
     1625  public:
     1626
     1627    /// \brief Constructor
     1628    ///
     1629    /// Constructor.
     1630    MaxWeightedMatching(const UGraph& ugraph, const WeightMap& weight)
     1631      : _ugraph(ugraph), _weight(weight), _matching(0),
     1632        _node_potential(0), _blossom_potential(), _blossom_node_list(),
     1633        _node_num(0), _blossom_num(0),
     1634
     1635        _blossom_index(0), _blossom_set(0), _blossom_data(0),
     1636        _node_index(0), _node_heap_index(0), _node_data(0),
     1637        _tree_set_index(0), _tree_set(0),
     1638
     1639        _delta1_index(0), _delta1(0),
     1640        _delta2_index(0), _delta2(0),
     1641        _delta3_index(0), _delta3(0),
     1642        _delta4_index(0), _delta4(0),
     1643
     1644        _delta_sum() {}
     1645
     1646    ~MaxWeightedMatching() {
     1647      destroyStructures();
     1648    }
     1649
     1650    /// \name Execution control
     1651    /// The simplest way to execute the algorithm is to use the member
     1652    /// \c run() member function.
     1653
     1654    ///@{
     1655
     1656    /// \brief Initialize the algorithm
     1657    ///
     1658    /// Initialize the algorithm
     1659    void init() {
     1660      createStructures();
     1661
     1662      for (EdgeIt e(_ugraph); e != INVALID; ++e) {
     1663        _node_heap_index->set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP);
     1664      }
     1665      for (NodeIt n(_ugraph); n != INVALID; ++n) {
     1666        _delta1_index->set(n, _delta1->PRE_HEAP);
     1667      }
     1668      for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
     1669        _delta3_index->set(e, _delta3->PRE_HEAP);
     1670      }
     1671      for (int i = 0; i < _blossom_num; ++i) {
     1672        _delta2_index->set(i, _delta2->PRE_HEAP);
     1673        _delta4_index->set(i, _delta4->PRE_HEAP);
     1674      }
     1675
     1676      int index = 0;
     1677      for (NodeIt n(_ugraph); n != INVALID; ++n) {
     1678        Value max = 0;
     1679        for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
     1680          if (_ugraph.target(e) == n) continue;
     1681          if ((dualScale * _weight[e]) / 2 > max) {
     1682            max = (dualScale * _weight[e]) / 2;
     1683          }
     1684        }
     1685        _node_index->set(n, index);
     1686        (*_node_data)[index].pot = max;
     1687        _delta1->push(n, max);
     1688        int blossom =
     1689          _blossom_set->insert(n, std::numeric_limits<Value>::max());
     1690
     1691        _tree_set->insert(blossom);
     1692
     1693        (*_blossom_data)[blossom].status = EVEN;
     1694        (*_blossom_data)[blossom].pred = INVALID;
     1695        (*_blossom_data)[blossom].next = INVALID;
     1696        (*_blossom_data)[blossom].pot = 0;
     1697        (*_blossom_data)[blossom].offset = 0;
     1698        ++index;
     1699      }
     1700      for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
     1701        int si = (*_node_index)[_ugraph.source(e)];
     1702        int ti = (*_node_index)[_ugraph.target(e)];
     1703        if (_ugraph.source(e) != _ugraph.target(e)) {
     1704          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
     1705                            dualScale * _weight[e]) / 2);
     1706        }
     1707      }
     1708    }
     1709
     1710    /// \brief Starts the algorithm
     1711    ///
     1712    /// Starts the algorithm
     1713    void start() {
     1714      enum OpType {
     1715        D1, D2, D3, D4
     1716      };
     1717
     1718      int unmatched = _node_num;
     1719      while (unmatched > 0) {
     1720        Value d1 = !_delta1->empty() ?
     1721          _delta1->prio() : std::numeric_limits<Value>::max();
     1722
     1723        Value d2 = !_delta2->empty() ?
     1724          _delta2->prio() : std::numeric_limits<Value>::max();
     1725
     1726        Value d3 = !_delta3->empty() ?
     1727          _delta3->prio() : std::numeric_limits<Value>::max();
     1728
     1729        Value d4 = !_delta4->empty() ?
     1730          _delta4->prio() : std::numeric_limits<Value>::max();
     1731
     1732        _delta_sum = d1; OpType ot = D1;
     1733        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
     1734        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
     1735        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
     1736
     1737       
     1738        switch (ot) {
     1739        case D1:
     1740          {
     1741            Node n = _delta1->top();
     1742            unmatchNode(n);
     1743            --unmatched;
     1744          }
     1745          break;
     1746        case D2:
     1747          {
     1748            int blossom = _delta2->top();
     1749            Node n = _blossom_set->classTop(blossom);
     1750            Edge e = (*_node_data)[(*_node_index)[n]].heap.top();
     1751            extendOnEdge(e);
     1752          }
     1753          break;
     1754        case D3:
     1755          {
     1756            UEdge e = _delta3->top();
     1757           
     1758            int left_blossom = _blossom_set->find(_ugraph.source(e));
     1759            int right_blossom = _blossom_set->find(_ugraph.target(e));
     1760         
     1761            if (left_blossom == right_blossom) {
     1762              _delta3->pop();
     1763            } else {
     1764              int left_tree;
     1765              if ((*_blossom_data)[left_blossom].status == EVEN) {
     1766                left_tree = _tree_set->find(left_blossom);
     1767              } else {
     1768                left_tree = -1;
     1769                ++unmatched;
     1770              }
     1771              int right_tree;
     1772              if ((*_blossom_data)[right_blossom].status == EVEN) {
     1773                right_tree = _tree_set->find(right_blossom);
     1774              } else {
     1775                right_tree = -1;
     1776                ++unmatched;
     1777              }
     1778           
     1779              if (left_tree == right_tree) {
     1780                shrinkOnEdge(e, left_tree);
     1781              } else {
     1782                augmentOnEdge(e);
     1783                unmatched -= 2;
     1784              }
     1785            }
     1786          } break;
     1787        case D4:
     1788          splitBlossom(_delta4->top());
     1789          break;
     1790        }
     1791      }
     1792      extractMatching();
     1793    }
     1794
     1795    /// \brief Runs %MaxWeightedMatching algorithm.
     1796    ///
     1797    /// This method runs the %MaxWeightedMatching algorithm.
     1798    ///
     1799    /// \note mwm.run() is just a shortcut of the following code.
     1800    /// \code
     1801    ///   mwm.init();
     1802    ///   mwm.start();
     1803    /// \endcode
     1804    void run() {
     1805      init();
     1806      start();
     1807    }
     1808
     1809    /// @}
     1810
     1811    /// \name Primal solution
     1812    /// Functions for get the primal solution, ie. the matching.
     1813   
     1814    /// @{
     1815
     1816    /// \brief Returns the matching value.
     1817    ///
     1818    /// Returns the matching value.
     1819    Value matchingValue() const {
     1820      Value sum = 0;
     1821      for (NodeIt n(_ugraph); n != INVALID; ++n) {
     1822        if ((*_matching)[n] != INVALID) {
     1823          sum += _weight[(*_matching)[n]];
     1824        }
     1825      }
     1826      return sum /= 2;
     1827    }
     1828
     1829    /// \brief Returns true when the edge is in the matching.
     1830    ///
     1831    /// Returns true when the edge is in the matching.
     1832    bool matching(const UEdge& edge) const {
     1833      return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true);
     1834    }
     1835
     1836    /// \brief Returns the incident matching edge.
     1837    ///
     1838    /// Returns the incident matching edge from given node. If the
     1839    /// node is not matched then it gives back \c INVALID.
     1840    Edge matching(const Node& node) const {
     1841      return (*_matching)[node];
     1842    }
     1843
     1844    /// \brief Returns the mate of the node.
     1845    ///
     1846    /// Returns the adjancent node in a mathcing edge. If the node is
     1847    /// not matched then it gives back \c INVALID.
     1848    Node mate(const Node& node) const {
     1849      return (*_matching)[node] != INVALID ?
     1850        _ugraph.target((*_matching)[node]) : INVALID;
     1851    }
     1852
     1853    /// @}
     1854
     1855    /// \name Dual solution
     1856    /// Functions for get the dual solution.
     1857   
     1858    /// @{
     1859
     1860    /// \brief Returns the value of the dual solution.
     1861    ///
     1862    /// Returns the value of the dual solution. It should be equal to
     1863    /// the primal value scaled by \ref dualScale "dual scale".
     1864    Value dualValue() const {
     1865      Value sum = 0;
     1866      for (NodeIt n(_ugraph); n != INVALID; ++n) {
     1867        sum += nodeValue(n);
     1868      }
     1869      for (int i = 0; i < blossomNum(); ++i) {
     1870        sum += blossomValue(i) * (blossomSize(i) / 2);
     1871      }
     1872      return sum;
     1873    }
     1874
     1875    /// \brief Returns the value of the node.
     1876    ///
     1877    /// Returns the the value of the node.
     1878    Value nodeValue(const Node& n) const {
     1879      return (*_node_potential)[n];
     1880    }
     1881
     1882    /// \brief Returns the number of the blossoms in the basis.
     1883    ///
     1884    /// Returns the number of the blossoms in the basis.
     1885    /// \see BlossomIt
     1886    int blossomNum() const {
     1887      return _blossom_potential.size();
     1888    }
     1889
     1890
     1891    /// \brief Returns the number of the nodes in the blossom.
     1892    ///
     1893    /// Returns the number of the nodes in the blossom.
     1894    int blossomSize(int k) const {
     1895      return _blossom_potential[k].end - _blossom_potential[k].begin;
     1896    }
     1897
     1898    /// \brief Returns the value of the blossom.
     1899    ///
     1900    /// Returns the the value of the blossom.
     1901    /// \see BlossomIt
     1902    Value blossomValue(int k) const {
     1903      return _blossom_potential[k].value;
     1904    }
     1905
     1906    /// \brief Lemon iterator for get the items of the blossom.
     1907    ///
     1908    /// Lemon iterator for get the nodes of the blossom. This class
     1909    /// provides a common style lemon iterator which gives back a
     1910    /// subset of the nodes.
     1911    class BlossomIt {
     1912    public:
     1913
     1914      /// \brief Constructor.
     1915      ///
     1916      /// Constructor for get the nodes of the variable.
     1917      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
     1918        : _algorithm(&algorithm)
     1919      {
     1920        _index = _algorithm->_blossom_potential[variable].begin;
     1921        _last = _algorithm->_blossom_potential[variable].end;
     1922      }
     1923
     1924      /// \brief Invalid constructor.
     1925      ///
     1926      /// Invalid constructor.
     1927      BlossomIt(Invalid) : _index(-1) {}
     1928
     1929      /// \brief Conversion to node.
     1930      ///
     1931      /// Conversion to node.
     1932      operator Node() const {
     1933        return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
     1934      }
     1935
     1936      /// \brief Increment operator.
     1937      ///
     1938      /// Increment operator.
     1939      BlossomIt& operator++() {
     1940        ++_index;
     1941        if (_index == _last) {
     1942          _index = -1;
     1943        }
     1944        return *this;
     1945      }
     1946
     1947      bool operator==(const BlossomIt& it) const {
     1948        return _index == it._index;
     1949      }
     1950      bool operator!=(const BlossomIt& it) const {
     1951        return _index != it._index;
     1952      }
     1953
     1954    private:
     1955      const MaxWeightedMatching* _algorithm;
     1956      int _last;
     1957      int _index;
     1958    };
     1959
     1960    /// @}
     1961   
     1962  };
     1963
     1964  /// \ingroup matching
     1965  ///
     1966  /// \brief Weighted perfect matching in general undirected graphs
     1967  ///
     1968  /// This class provides an efficient implementation of Edmond's
     1969  /// maximum weighted perfecr matching algorithm. The implementation
     1970  /// is based on extensive use of priority queues and provides
     1971  /// \f$O(nm\log(n))\f$ time complexity.
     1972  ///
     1973  /// The maximum weighted matching problem is to find undirected
     1974  /// edges in the graph with maximum overall weight and no two of
     1975  /// them shares their endpoints and covers all nodes. The problem
     1976  /// can be formulated with the next linear program:
     1977  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
     1978  ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
     1979  /// \f[x_e \ge 0\quad \forall e\in E\f]
     1980  /// \f[\max \sum_{e\in E}x_ew_e\f]
     1981  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
     1982  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in
     1983  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
     1984  /// the nodes.
     1985  ///
     1986  /// The algorithm calculates an optimal matching and a proof of the
     1987  /// optimality. The solution of the dual problem can be used to check
     1988  /// the result of the algorithm. The dual linear problem is the next:
     1989  /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
     1990  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
     1991  /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
     1992  ///
     1993  /// The algorithm can be executed with \c run() or the \c init() and
     1994  /// then the \c start() member functions. After it the matching can
     1995  /// be asked with \c matching() or mate() functions. The dual
     1996  /// solution can be get with \c nodeValue(), \c blossomNum() and \c
     1997  /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
     1998  /// "BlossomIt" nested class which is able to iterate on the nodes
     1999  /// of a blossom. If the value type is integral then the dual
     2000  /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
     2001  ///
     2002  /// \author Balazs Dezso
     2003  template <typename _UGraph,
     2004            typename _WeightMap = typename _UGraph::template UEdgeMap<int> >
     2005  class MaxWeightedPerfectMatching {
     2006  public:
     2007
     2008    typedef _UGraph UGraph;
     2009    typedef _WeightMap WeightMap;
     2010    typedef typename WeightMap::Value Value;
     2011
     2012    /// \brief Scaling factor for dual solution
     2013    ///
     2014    /// Scaling factor for dual solution, it is equal to 4 or 1
     2015    /// according to the value type.
     2016    static const int dualScale =
     2017      std::numeric_limits<Value>::is_integer ? 4 : 1;
     2018
     2019    typedef typename UGraph::template NodeMap<typename UGraph::Edge>
     2020    MatchingMap;
     2021   
     2022  private:
     2023
     2024    UGRAPH_TYPEDEFS(typename UGraph);
     2025
     2026    typedef typename UGraph::template NodeMap<Value> NodePotential;
     2027    typedef std::vector<Node> BlossomNodeList;
     2028
     2029    struct BlossomVariable {
     2030      int begin, end;
     2031      Value value;
     2032     
     2033      BlossomVariable(int _begin, int _end, Value _value)
     2034        : begin(_begin), end(_end), value(_value) {}
     2035
     2036    };
     2037
     2038    typedef std::vector<BlossomVariable> BlossomPotential;
     2039
     2040    const UGraph& _ugraph;
     2041    const WeightMap& _weight;
     2042
     2043    MatchingMap* _matching;
     2044
     2045    NodePotential* _node_potential;
     2046
     2047    BlossomPotential _blossom_potential;
     2048    BlossomNodeList _blossom_node_list;
     2049
     2050    int _node_num;
     2051    int _blossom_num;
     2052
     2053    typedef typename UGraph::template NodeMap<int> NodeIntMap;
     2054    typedef typename UGraph::template EdgeMap<int> EdgeIntMap;
     2055    typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap;
     2056    typedef IntegerMap<int> IntIntMap;
     2057
     2058    enum Status {
     2059      EVEN = -1, MATCHED = 0, ODD = 1
     2060    };
     2061
     2062    typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
     2063    struct BlossomData {
     2064      int tree;
     2065      Status status;
     2066      Edge pred, next;
     2067      Value pot, offset;
     2068    };
     2069
     2070    NodeIntMap *_blossom_index;
     2071    BlossomSet *_blossom_set;
     2072    IntegerMap<BlossomData>* _blossom_data;
     2073
     2074    NodeIntMap *_node_index;
     2075    EdgeIntMap *_node_heap_index;
     2076
     2077    struct NodeData {
     2078     
     2079      NodeData(EdgeIntMap& node_heap_index)
     2080        : heap(node_heap_index) {}
     2081     
     2082      int blossom;
     2083      Value pot;
     2084      BinHeap<Value, EdgeIntMap> heap;
     2085      std::map<int, Edge> heap_index;
     2086     
     2087      int tree;
     2088    };
     2089
     2090    IntegerMap<NodeData>* _node_data;
     2091
     2092    typedef ExtendFindEnum<IntIntMap> TreeSet;
     2093
     2094    IntIntMap *_tree_set_index;
     2095    TreeSet *_tree_set;
     2096
     2097    IntIntMap *_delta2_index;
     2098    BinHeap<Value, IntIntMap> *_delta2;
     2099   
     2100    UEdgeIntMap *_delta3_index;
     2101    BinHeap<Value, UEdgeIntMap> *_delta3;
     2102
     2103    IntIntMap *_delta4_index;
     2104    BinHeap<Value, IntIntMap> *_delta4;
     2105
     2106    Value _delta_sum;
     2107
     2108    void createStructures() {
     2109      _node_num = countNodes(_ugraph);
     2110      _blossom_num = _node_num * 3 / 2;
     2111
     2112      if (!_matching) {
     2113        _matching = new MatchingMap(_ugraph);
     2114      }
     2115      if (!_node_potential) {
     2116        _node_potential = new NodePotential(_ugraph);
     2117      }
     2118      if (!_blossom_set) {
     2119        _blossom_index = new NodeIntMap(_ugraph);
     2120        _blossom_set = new BlossomSet(*_blossom_index);
     2121        _blossom_data = new IntegerMap<BlossomData>(_blossom_num);
     2122      }
     2123
     2124      if (!_node_index) {
     2125        _node_index = new NodeIntMap(_ugraph);
     2126        _node_heap_index = new EdgeIntMap(_ugraph);
     2127        _node_data = new IntegerMap<NodeData>(_node_num,
     2128                                              NodeData(*_node_heap_index));
     2129      }
     2130
     2131      if (!_tree_set) {
     2132        _tree_set_index = new IntIntMap(_blossom_num);
     2133        _tree_set = new TreeSet(*_tree_set_index);
     2134      }
     2135      if (!_delta2) {
     2136        _delta2_index = new IntIntMap(_blossom_num);
     2137        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
     2138      }
     2139      if (!_delta3) {
     2140        _delta3_index = new UEdgeIntMap(_ugraph);
     2141        _delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index);
     2142      }
     2143      if (!_delta4) {
     2144        _delta4_index = new IntIntMap(_blossom_num);
     2145        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
     2146      }
     2147    }
     2148
     2149    void destroyStructures() {
     2150      _node_num = countNodes(_ugraph);
     2151      _blossom_num = _node_num * 3 / 2;
     2152
     2153      if (_matching) {
     2154        delete _matching;
     2155      }
     2156      if (_node_potential) {
     2157        delete _node_potential;
     2158      }
     2159      if (_blossom_set) {
     2160        delete _blossom_index;
     2161        delete _blossom_set;
     2162        delete _blossom_data;
     2163      }
     2164
     2165      if (_node_index) {
     2166        delete _node_index;
     2167        delete _node_heap_index;
     2168        delete _node_data;                           
     2169      }
     2170
     2171      if (_tree_set) {
     2172        delete _tree_set_index;
     2173        delete _tree_set;
     2174      }
     2175      if (_delta2) {
     2176        delete _delta2_index;
     2177        delete _delta2;
     2178      }
     2179      if (_delta3) {
     2180        delete _delta3_index;
     2181        delete _delta3;
     2182      }
     2183      if (_delta4) {
     2184        delete _delta4_index;
     2185        delete _delta4;
     2186      }
     2187    }
     2188
     2189    void matchedToEven(int blossom, int tree) {
     2190      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
     2191        _delta2->erase(blossom);
     2192      }
     2193
     2194      if (!_blossom_set->trivial(blossom)) {
     2195        (*_blossom_data)[blossom].pot -=
     2196          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
     2197      }
     2198
     2199      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
     2200           n != INVALID; ++n) {
     2201
     2202        _blossom_set->increase(n, std::numeric_limits<Value>::max());
     2203        int ni = (*_node_index)[n];
     2204
     2205        (*_node_data)[ni].heap.clear();
     2206        (*_node_data)[ni].heap_index.clear();
     2207
     2208        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
     2209
     2210        for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
     2211          Node v = _ugraph.source(e);
     2212          int vb = _blossom_set->find(v);
     2213          int vi = (*_node_index)[v];
     2214
     2215          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
     2216            dualScale * _weight[e];
     2217
     2218          if ((*_blossom_data)[vb].status == EVEN) {
     2219            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
     2220              _delta3->push(e, rw / 2);
     2221            }
     2222          } else {
     2223            typename std::map<int, Edge>::iterator it =
     2224              (*_node_data)[vi].heap_index.find(tree);   
     2225
     2226            if (it != (*_node_data)[vi].heap_index.end()) {
     2227              if ((*_node_data)[vi].heap[it->second] > rw) {
     2228                (*_node_data)[vi].heap.replace(it->second, e);
     2229                (*_node_data)[vi].heap.decrease(e, rw);
     2230                it->second = e;
     2231              }
     2232            } else {
     2233              (*_node_data)[vi].heap.push(e, rw);
     2234              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
     2235            }
     2236
     2237            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
     2238              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
     2239
     2240              if ((*_blossom_data)[vb].status == MATCHED) {
     2241                if (_delta2->state(vb) != _delta2->IN_HEAP) {
     2242                  _delta2->push(vb, _blossom_set->classPrio(vb) -
     2243                               (*_blossom_data)[vb].offset);
     2244                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
     2245                           (*_blossom_data)[vb].offset){
     2246                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
     2247                                   (*_blossom_data)[vb].offset);
     2248                }
     2249              }
     2250            }
     2251          }
     2252        }
     2253      }
     2254      (*_blossom_data)[blossom].offset = 0;
     2255    }
     2256
     2257    void matchedToOdd(int blossom) {
     2258      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
     2259        _delta2->erase(blossom);
     2260      }
     2261      (*_blossom_data)[blossom].offset += _delta_sum;
     2262      if (!_blossom_set->trivial(blossom)) {
     2263        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
     2264                     (*_blossom_data)[blossom].offset);
     2265      }
     2266    }
     2267
     2268    void evenToMatched(int blossom, int tree) {
     2269      if (!_blossom_set->trivial(blossom)) {
     2270        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
     2271      }
     2272
     2273      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
     2274           n != INVALID; ++n) {
     2275        int ni = (*_node_index)[n];
     2276        (*_node_data)[ni].pot -= _delta_sum;
     2277
     2278        for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
     2279          Node v = _ugraph.source(e);
     2280          int vb = _blossom_set->find(v);
     2281          int vi = (*_node_index)[v];
     2282
     2283          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
     2284            dualScale * _weight[e];
     2285
     2286          if (vb == blossom) {
     2287            if (_delta3->state(e) == _delta3->IN_HEAP) {
     2288              _delta3->erase(e);
     2289            }
     2290          } else if ((*_blossom_data)[vb].status == EVEN) {
     2291
     2292            if (_delta3->state(e) == _delta3->IN_HEAP) {
     2293              _delta3->erase(e);
     2294            }
     2295
     2296            int vt = _tree_set->find(vb);
     2297           
     2298            if (vt != tree) {
     2299
     2300              Edge r = _ugraph.oppositeEdge(e);
     2301
     2302              typename std::map<int, Edge>::iterator it =
     2303                (*_node_data)[ni].heap_index.find(vt);   
     2304
     2305              if (it != (*_node_data)[ni].heap_index.end()) {
     2306                if ((*_node_data)[ni].heap[it->second] > rw) {
     2307                  (*_node_data)[ni].heap.replace(it->second, r);
     2308                  (*_node_data)[ni].heap.decrease(r, rw);
     2309                  it->second = r;
     2310                }
     2311              } else {
     2312                (*_node_data)[ni].heap.push(r, rw);
     2313                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
     2314              }
     2315             
     2316              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
     2317                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
     2318               
     2319                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
     2320                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
     2321                               (*_blossom_data)[blossom].offset);
     2322                } else if ((*_delta2)[blossom] >
     2323                           _blossom_set->classPrio(blossom) -
     2324                           (*_blossom_data)[blossom].offset){
     2325                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
     2326                                   (*_blossom_data)[blossom].offset);
     2327                }
     2328              }
     2329            }
     2330          } else {
     2331         
     2332            typename std::map<int, Edge>::iterator it =
     2333              (*_node_data)[vi].heap_index.find(tree);
     2334
     2335            if (it != (*_node_data)[vi].heap_index.end()) {
     2336              (*_node_data)[vi].heap.erase(it->second);
     2337              (*_node_data)[vi].heap_index.erase(it);
     2338              if ((*_node_data)[vi].heap.empty()) {
     2339                _blossom_set->increase(v, std::numeric_limits<Value>::max());
     2340              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
     2341                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
     2342              }
     2343             
     2344              if ((*_blossom_data)[vb].status == MATCHED) {
     2345                if (_blossom_set->classPrio(vb) ==
     2346                    std::numeric_limits<Value>::max()) {
     2347                  _delta2->erase(vb);
     2348                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
     2349                           (*_blossom_data)[vb].offset) {
     2350                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
     2351                                   (*_blossom_data)[vb].offset);
     2352                }
     2353              }
     2354            }   
     2355          }
     2356        }
     2357      }
     2358    }
     2359
     2360    void oddToMatched(int blossom) {
     2361      (*_blossom_data)[blossom].offset -= _delta_sum;
     2362
     2363      if (_blossom_set->classPrio(blossom) !=
     2364          std::numeric_limits<Value>::max()) {
     2365        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
     2366                       (*_blossom_data)[blossom].offset);
     2367      }
     2368
     2369      if (!_blossom_set->trivial(blossom)) {
     2370        _delta4->erase(blossom);
     2371      }
     2372    }
     2373
     2374    void oddToEven(int blossom, int tree) {
     2375      if (!_blossom_set->trivial(blossom)) {
     2376        _delta4->erase(blossom);
     2377        (*_blossom_data)[blossom].pot -=
     2378          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
     2379      }
     2380
     2381      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
     2382           n != INVALID; ++n) {
     2383        int ni = (*_node_index)[n];
     2384
     2385        _blossom_set->increase(n, std::numeric_limits<Value>::max());
     2386
     2387        (*_node_data)[ni].heap.clear();
     2388        (*_node_data)[ni].heap_index.clear();
     2389        (*_node_data)[ni].pot +=
     2390          2 * _delta_sum - (*_blossom_data)[blossom].offset;
     2391
     2392        for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
     2393          Node v = _ugraph.source(e);
     2394          int vb = _blossom_set->find(v);
     2395          int vi = (*_node_index)[v];
     2396
     2397          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
     2398            dualScale * _weight[e];
     2399
     2400          if ((*_blossom_data)[vb].status == EVEN) {
     2401            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
     2402              _delta3->push(e, rw / 2);
     2403            }
     2404          } else {
     2405         
     2406            typename std::map<int, Edge>::iterator it =
     2407              (*_node_data)[vi].heap_index.find(tree);   
     2408
     2409            if (it != (*_node_data)[vi].heap_index.end()) {
     2410              if ((*_node_data)[vi].heap[it->second] > rw) {
     2411                (*_node_data)[vi].heap.replace(it->second, e);
     2412                (*_node_data)[vi].heap.decrease(e, rw);
     2413                it->second = e;
     2414              }
     2415            } else {
     2416              (*_node_data)[vi].heap.push(e, rw);
     2417              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
     2418            }
     2419
     2420            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
     2421              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
     2422
     2423              if ((*_blossom_data)[vb].status == MATCHED) {
     2424                if (_delta2->state(vb) != _delta2->IN_HEAP) {
     2425                  _delta2->push(vb, _blossom_set->classPrio(vb) -
     2426                               (*_blossom_data)[vb].offset);
     2427                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
     2428                           (*_blossom_data)[vb].offset) {
     2429                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
     2430                                   (*_blossom_data)[vb].offset);
     2431                }
     2432              }
     2433            }
     2434          }
     2435        }
     2436      }
     2437      (*_blossom_data)[blossom].offset = 0;
     2438    }
     2439   
     2440    void alternatePath(int even, int tree) {
     2441      int odd;
     2442
     2443      evenToMatched(even, tree);
     2444      (*_blossom_data)[even].status = MATCHED;
     2445
     2446      while ((*_blossom_data)[even].pred != INVALID) {
     2447        odd = _blossom_set->find(_ugraph.target((*_blossom_data)[even].pred));
     2448        (*_blossom_data)[odd].status = MATCHED;
     2449        oddToMatched(odd);
     2450        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
     2451     
     2452        even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred));
     2453        (*_blossom_data)[even].status = MATCHED;
     2454        evenToMatched(even, tree);
     2455        (*_blossom_data)[even].next =
     2456          _ugraph.oppositeEdge((*_blossom_data)[odd].pred);
     2457      }
     2458     
     2459    }
     2460
     2461    void destroyTree(int tree) {
     2462      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
     2463        if ((*_blossom_data)[b].status == EVEN) {
     2464          (*_blossom_data)[b].status = MATCHED;
     2465          evenToMatched(b, tree);
     2466        } else if ((*_blossom_data)[b].status == ODD) {
     2467          (*_blossom_data)[b].status = MATCHED;
     2468          oddToMatched(b);
     2469        }
     2470      }
     2471      _tree_set->eraseClass(tree);
     2472    }
     2473
     2474    void augmentOnEdge(const UEdge& edge) {
     2475     
     2476      int left = _blossom_set->find(_ugraph.source(edge));
     2477      int right = _blossom_set->find(_ugraph.target(edge));
     2478
     2479      int left_tree = _tree_set->find(left);
     2480      alternatePath(left, left_tree);
     2481      destroyTree(left_tree);
     2482
     2483      int right_tree = _tree_set->find(right);
     2484      alternatePath(right, right_tree);
     2485      destroyTree(right_tree);
     2486
     2487      (*_blossom_data)[left].next = _ugraph.direct(edge, true);
     2488      (*_blossom_data)[right].next = _ugraph.direct(edge, false);
     2489    }
     2490
     2491    void extendOnEdge(const Edge& edge) {
     2492      int base = _blossom_set->find(_ugraph.target(edge));
     2493      int tree = _tree_set->find(base);
     2494     
     2495      int odd = _blossom_set->find(_ugraph.source(edge));
     2496      _tree_set->insert(odd, tree);
     2497      (*_blossom_data)[odd].status = ODD;
     2498      matchedToOdd(odd);
     2499      (*_blossom_data)[odd].pred = edge;
     2500
     2501      int even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].next));
     2502      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
     2503      _tree_set->insert(even, tree);
     2504      (*_blossom_data)[even].status = EVEN;
     2505      matchedToEven(even, tree);
     2506    }
     2507   
     2508    void shrinkOnEdge(const UEdge& uedge, int tree) {
     2509      int nca = -1;
     2510      std::vector<int> left_path, right_path;
     2511
     2512      {
     2513        std::set<int> left_set, right_set;
     2514        int left = _blossom_set->find(_ugraph.source(uedge));
     2515        left_path.push_back(left);
     2516        left_set.insert(left);
     2517
     2518        int right = _blossom_set->find(_ugraph.target(uedge));
     2519        right_path.push_back(right);
     2520        right_set.insert(right);
     2521
     2522        while (true) {
     2523
     2524          if ((*_blossom_data)[left].pred == INVALID) break;
     2525
     2526          left =
     2527            _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
     2528          left_path.push_back(left);
     2529          left =
     2530            _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
     2531          left_path.push_back(left);
     2532
     2533          left_set.insert(left);
     2534
     2535          if (right_set.find(left) != right_set.end()) {
     2536            nca = left;
     2537            break;
     2538          }
     2539
     2540          if ((*_blossom_data)[right].pred == INVALID) break;
     2541
     2542          right =
     2543            _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
     2544          right_path.push_back(right);
     2545          right =
     2546            _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
     2547          right_path.push_back(right);
     2548
     2549          right_set.insert(right);
     2550 
     2551          if (left_set.find(right) != left_set.end()) {
     2552            nca = right;
     2553            break;
     2554          }
     2555
     2556        }
     2557
     2558        if (nca == -1) {
     2559          if ((*_blossom_data)[left].pred == INVALID) {
     2560            nca = right;
     2561            while (left_set.find(nca) == left_set.end()) {
     2562              nca =
     2563                _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
     2564              right_path.push_back(nca);
     2565              nca =
     2566                _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
     2567              right_path.push_back(nca);
     2568            }
     2569          } else {
     2570            nca = left;
     2571            while (right_set.find(nca) == right_set.end()) {
     2572              nca =
     2573                _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
     2574              left_path.push_back(nca);
     2575              nca =
     2576                _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
     2577              left_path.push_back(nca);
     2578            }
     2579          }
     2580        }
     2581      }
     2582
     2583      std::vector<int> subblossoms;
     2584      Edge prev;
     2585
     2586      prev = _ugraph.direct(uedge, true);
     2587      for (int i = 0; left_path[i] != nca; i += 2) {
     2588        subblossoms.push_back(left_path[i]);
     2589        (*_blossom_data)[left_path[i]].next = prev;
     2590        _tree_set->erase(left_path[i]);
     2591
     2592        subblossoms.push_back(left_path[i + 1]);
     2593        (*_blossom_data)[left_path[i + 1]].status = EVEN;
     2594        oddToEven(left_path[i + 1], tree);
     2595        _tree_set->erase(left_path[i + 1]);
     2596        prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred);
     2597      }
     2598
     2599      int k = 0;
     2600      while (right_path[k] != nca) ++k;
     2601
     2602      subblossoms.push_back(nca);
     2603      (*_blossom_data)[nca].next = prev;
     2604     
     2605      for (int i = k - 2; i >= 0; i -= 2) {
     2606        subblossoms.push_back(right_path[i + 1]);
     2607        (*_blossom_data)[right_path[i + 1]].status = EVEN;
     2608        oddToEven(right_path[i + 1], tree);
     2609        _tree_set->erase(right_path[i + 1]);
     2610
     2611        (*_blossom_data)[right_path[i + 1]].next =
     2612          (*_blossom_data)[right_path[i + 1]].pred;
     2613
     2614        subblossoms.push_back(right_path[i]);
     2615        _tree_set->erase(right_path[i]);
     2616      }
     2617
     2618      int surface =
     2619        _blossom_set->join(subblossoms.begin(), subblossoms.end());
     2620
     2621      for (int i = 0; i < int(subblossoms.size()); ++i) {
     2622        if (!_blossom_set->trivial(subblossoms[i])) {
     2623          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
     2624        }
     2625        (*_blossom_data)[subblossoms[i]].status = MATCHED;
     2626      }
     2627
     2628      (*_blossom_data)[surface].pot = -2 * _delta_sum;
     2629      (*_blossom_data)[surface].offset = 0;
     2630      (*_blossom_data)[surface].status = EVEN;
     2631      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
     2632      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
     2633
     2634      _tree_set->insert(surface, tree);
     2635      _tree_set->erase(nca);
     2636    }
     2637
     2638    void splitBlossom(int blossom) {
     2639      Edge next = (*_blossom_data)[blossom].next;
     2640      Edge pred = (*_blossom_data)[blossom].pred;
     2641
     2642      int tree = _tree_set->find(blossom);
     2643
     2644      (*_blossom_data)[blossom].status = MATCHED;
     2645      oddToMatched(blossom);
     2646      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
     2647        _delta2->erase(blossom);
     2648      }
     2649
     2650      std::vector<int> subblossoms;
     2651      _blossom_set->split(blossom, std::back_inserter(subblossoms));
     2652
     2653      Value offset = (*_blossom_data)[blossom].offset;
     2654      int b = _blossom_set->find(_ugraph.source(pred));
     2655      int d = _blossom_set->find(_ugraph.source(next));
     2656     
     2657      int ib, id;
     2658      for (int i = 0; i < int(subblossoms.size()); ++i) {
     2659        if (subblossoms[i] == b) ib = i;
     2660        if (subblossoms[i] == d) id = i;
     2661
     2662        (*_blossom_data)[subblossoms[i]].offset = offset;
     2663        if (!_blossom_set->trivial(subblossoms[i])) {
     2664          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
     2665        }
     2666        if (_blossom_set->classPrio(subblossoms[i]) !=
     2667            std::numeric_limits<Value>::max()) {
     2668          _delta2->push(subblossoms[i],
     2669                        _blossom_set->classPrio(subblossoms[i]) -
     2670                        (*_blossom_data)[subblossoms[i]].offset);
     2671        }
     2672      }
     2673     
     2674      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
     2675        for (int i = (id + 1) % subblossoms.size();
     2676             i != ib; i = (i + 2) % subblossoms.size()) {
     2677          int sb = subblossoms[i];
     2678          int tb = subblossoms[(i + 1) % subblossoms.size()];
     2679          (*_blossom_data)[sb].next =
     2680            _ugraph.oppositeEdge((*_blossom_data)[tb].next);
     2681        }
     2682
     2683        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
     2684          int sb = subblossoms[i];
     2685          int tb = subblossoms[(i + 1) % subblossoms.size()];
     2686          int ub = subblossoms[(i + 2) % subblossoms.size()];
     2687
     2688          (*_blossom_data)[sb].status = ODD;
     2689          matchedToOdd(sb);
     2690          _tree_set->insert(sb, tree);
     2691          (*_blossom_data)[sb].pred = pred;
     2692          (*_blossom_data)[sb].next =
     2693                           _ugraph.oppositeEdge((*_blossom_data)[tb].next);
     2694         
     2695          pred = (*_blossom_data)[ub].next;
     2696     
     2697          (*_blossom_data)[tb].status = EVEN;
     2698          matchedToEven(tb, tree);
     2699          _tree_set->insert(tb, tree);
     2700          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
     2701        }
     2702     
     2703        (*_blossom_data)[subblossoms[id]].status = ODD;
     2704        matchedToOdd(subblossoms[id]);
     2705        _tree_set->insert(subblossoms[id], tree);
     2706        (*_blossom_data)[subblossoms[id]].next = next;
     2707        (*_blossom_data)[subblossoms[id]].pred = pred;
     2708     
     2709      } else {
     2710
     2711        for (int i = (ib + 1) % subblossoms.size();
     2712             i != id; i = (i + 2) % subblossoms.size()) {
     2713          int sb = subblossoms[i];
     2714          int tb = subblossoms[(i + 1) % subblossoms.size()];
     2715          (*_blossom_data)[sb].next =
     2716            _ugraph.oppositeEdge((*_blossom_data)[tb].next);
     2717        }       
     2718
     2719        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
     2720          int sb = subblossoms[i];
     2721          int tb = subblossoms[(i + 1) % subblossoms.size()];
     2722          int ub = subblossoms[(i + 2) % subblossoms.size()];
     2723
     2724          (*_blossom_data)[sb].status = ODD;
     2725          matchedToOdd(sb);
     2726          _tree_set->insert(sb, tree);
     2727          (*_blossom_data)[sb].next = next;
     2728          (*_blossom_data)[sb].pred =
     2729            _ugraph.oppositeEdge((*_blossom_data)[tb].next);
     2730
     2731          (*_blossom_data)[tb].status = EVEN;
     2732          matchedToEven(tb, tree);
     2733          _tree_set->insert(tb, tree);
     2734          (*_blossom_data)[tb].pred =
     2735            (*_blossom_data)[tb].next =
     2736            _ugraph.oppositeEdge((*_blossom_data)[ub].next);
     2737          next = (*_blossom_data)[ub].next;
     2738        }
     2739
     2740        (*_blossom_data)[subblossoms[ib]].status = ODD;
     2741        matchedToOdd(subblossoms[ib]);
     2742        _tree_set->insert(subblossoms[ib], tree);
     2743        (*_blossom_data)[subblossoms[ib]].next = next;
     2744        (*_blossom_data)[subblossoms[ib]].pred = pred;
     2745      }
     2746      _tree_set->erase(blossom);
     2747    }
     2748
     2749    void extractBlossom(int blossom, const Node& base, const Edge& matching) {
     2750      if (_blossom_set->trivial(blossom)) {
     2751        int bi = (*_node_index)[base];
     2752        Value pot = (*_node_data)[bi].pot;
     2753
     2754        _matching->set(base, matching);
     2755        _blossom_node_list.push_back(base);
     2756        _node_potential->set(base, pot);
     2757      } else {
     2758
     2759        Value pot = (*_blossom_data)[blossom].pot;
     2760        int bn = _blossom_node_list.size();
     2761       
     2762        std::vector<int> subblossoms;
     2763        _blossom_set->split(blossom, std::back_inserter(subblossoms));
     2764        int b = _blossom_set->find(base);
     2765        int ib = -1;
     2766        for (int i = 0; i < int(subblossoms.size()); ++i) {
     2767          if (subblossoms[i] == b) { ib = i; break; }
     2768        }
     2769       
     2770        for (int i = 1; i < int(subblossoms.size()); i += 2) {
     2771          int sb = subblossoms[(ib + i) % subblossoms.size()];
     2772          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
     2773         
     2774          Edge m = (*_blossom_data)[tb].next;
     2775          extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m));
     2776          extractBlossom(tb, _ugraph.source(m), m);
     2777        }
     2778        extractBlossom(subblossoms[ib], base, matching);     
     2779       
     2780        int en = _blossom_node_list.size();
     2781       
     2782        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
     2783      }
     2784    }
     2785
     2786    void extractMatching() {
     2787      std::vector<int> blossoms;
     2788      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
     2789        blossoms.push_back(c);
     2790      }
     2791
     2792      for (int i = 0; i < int(blossoms.size()); ++i) {
     2793
     2794        Value offset = (*_blossom_data)[blossoms[i]].offset;
     2795        (*_blossom_data)[blossoms[i]].pot += 2 * offset;
     2796        for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
     2797             n != INVALID; ++n) {
     2798          (*_node_data)[(*_node_index)[n]].pot -= offset;
     2799        }
     2800       
     2801        Edge matching = (*_blossom_data)[blossoms[i]].next;
     2802        Node base = _ugraph.source(matching);
     2803        extractBlossom(blossoms[i], base, matching);
     2804      }
     2805    }
     2806   
     2807  public:
     2808
     2809    /// \brief Constructor
     2810    ///
     2811    /// Constructor.
     2812    MaxWeightedPerfectMatching(const UGraph& ugraph, const WeightMap& weight)
     2813      : _ugraph(ugraph), _weight(weight), _matching(0),
     2814        _node_potential(0), _blossom_potential(), _blossom_node_list(),
     2815        _node_num(0), _blossom_num(0),
     2816
     2817        _blossom_index(0), _blossom_set(0), _blossom_data(0),
     2818        _node_index(0), _node_heap_index(0), _node_data(0),
     2819        _tree_set_index(0), _tree_set(0),
     2820
     2821        _delta2_index(0), _delta2(0),
     2822        _delta3_index(0), _delta3(0),
     2823        _delta4_index(0), _delta4(0),
     2824
     2825        _delta_sum() {}
     2826
     2827    ~MaxWeightedPerfectMatching() {
     2828      destroyStructures();
     2829    }
     2830
     2831    /// \name Execution control
     2832    /// The simplest way to execute the algorithm is to use the member
     2833    /// \c run() member function.
     2834
     2835    ///@{
     2836
     2837    /// \brief Initialize the algorithm
     2838    ///
     2839    /// Initialize the algorithm
     2840    void init() {
     2841      createStructures();
     2842
     2843      for (EdgeIt e(_ugraph); e != INVALID; ++e) {
     2844        _node_heap_index->set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP);
     2845      }
     2846      for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
     2847        _delta3_index->set(e, _delta3->PRE_HEAP);
     2848      }
     2849      for (int i = 0; i < _blossom_num; ++i) {
     2850        _delta2_index->set(i, _delta2->PRE_HEAP);
     2851        _delta4_index->set(i, _delta4->PRE_HEAP);
     2852      }
     2853
     2854      int index = 0;
     2855      for (NodeIt n(_ugraph); n != INVALID; ++n) {
     2856        Value max = std::numeric_limits<Value>::min();
     2857        for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
     2858          if (_ugraph.target(e) == n) continue;
     2859          if ((dualScale * _weight[e]) / 2 > max) {
     2860            max = (dualScale * _weight[e]) / 2;
     2861          }
     2862        }
     2863        _node_index->set(n, index);
     2864        (*_node_data)[index].pot = max;
     2865        int blossom =
     2866          _blossom_set->insert(n, std::numeric_limits<Value>::max());
     2867
     2868        _tree_set->insert(blossom);
     2869
     2870        (*_blossom_data)[blossom].status = EVEN;
     2871        (*_blossom_data)[blossom].pred = INVALID;
     2872        (*_blossom_data)[blossom].next = INVALID;
     2873        (*_blossom_data)[blossom].pot = 0;
     2874        (*_blossom_data)[blossom].offset = 0;
     2875        ++index;
     2876      }
     2877      for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
     2878        int si = (*_node_index)[_ugraph.source(e)];
     2879        int ti = (*_node_index)[_ugraph.target(e)];
     2880        if (_ugraph.source(e) != _ugraph.target(e)) {
     2881          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
     2882                            dualScale * _weight[e]) / 2);
     2883        }
     2884      }
     2885    }
     2886
     2887    /// \brief Starts the algorithm
     2888    ///
     2889    /// Starts the algorithm
     2890    bool start() {
     2891      enum OpType {
     2892        D2, D3, D4
     2893      };
     2894
     2895      int unmatched = _node_num;
     2896      while (unmatched > 0) {
     2897        Value d2 = !_delta2->empty() ?
     2898          _delta2->prio() : std::numeric_limits<Value>::max();
     2899
     2900        Value d3 = !_delta3->empty() ?
     2901          _delta3->prio() : std::numeric_limits<Value>::max();
     2902
     2903        Value d4 = !_delta4->empty() ?
     2904          _delta4->prio() : std::numeric_limits<Value>::max();
     2905
     2906        _delta_sum = d2; OpType ot = D2;
     2907        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
     2908        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
     2909
     2910        if (_delta_sum == std::numeric_limits<Value>::max()) {
     2911          return false;
     2912        }
     2913       
     2914        switch (ot) {
     2915        case D2:
     2916          {
     2917            int blossom = _delta2->top();
     2918            Node n = _blossom_set->classTop(blossom);
     2919            Edge e = (*_node_data)[(*_node_index)[n]].heap.top();
     2920            extendOnEdge(e);
     2921          }
     2922          break;
     2923        case D3:
     2924          {
     2925            UEdge e = _delta3->top();
     2926           
     2927            int left_blossom = _blossom_set->find(_ugraph.source(e));
     2928            int right_blossom = _blossom_set->find(_ugraph.target(e));
     2929         
     2930            if (left_blossom == right_blossom) {
     2931              _delta3->pop();
     2932            } else {
     2933              int left_tree = _tree_set->find(left_blossom);
     2934              int right_tree = _tree_set->find(right_blossom);
     2935           
     2936              if (left_tree == right_tree) {
     2937                shrinkOnEdge(e, left_tree);
     2938              } else {
     2939                augmentOnEdge(e);
     2940                unmatched -= 2;
     2941              }
     2942            }
     2943          } break;
     2944        case D4:
     2945          splitBlossom(_delta4->top());
     2946          break;
     2947        }
     2948      }
     2949      extractMatching();
     2950      return true;
     2951    }
     2952
     2953    /// \brief Runs %MaxWeightedPerfectMatching algorithm.
     2954    ///
     2955    /// This method runs the %MaxWeightedPerfectMatching algorithm.
     2956    ///
     2957    /// \note mwm.run() is just a shortcut of the following code.
     2958    /// \code
     2959    ///   mwm.init();
     2960    ///   mwm.start();
     2961    /// \endcode
     2962    bool run() {
     2963      init();
     2964      return start();
     2965    }
     2966
     2967    /// @}
     2968
     2969    /// \name Primal solution
     2970    /// Functions for get the primal solution, ie. the matching.
     2971   
     2972    /// @{
     2973
     2974    /// \brief Returns the matching value.
     2975    ///
     2976    /// Returns the matching value.
     2977    Value matchingValue() const {
     2978      Value sum = 0;
     2979      for (NodeIt n(_ugraph); n != INVALID; ++n) {
     2980        if ((*_matching)[n] != INVALID) {
     2981          sum += _weight[(*_matching)[n]];
     2982        }
     2983      }
     2984      return sum /= 2;
     2985    }
     2986
     2987    /// \brief Returns true when the edge is in the matching.
     2988    ///
     2989    /// Returns true when the edge is in the matching.
     2990    bool matching(const UEdge& edge) const {
     2991      return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true);
     2992    }
     2993
     2994    /// \brief Returns the incident matching edge.
     2995    ///
     2996    /// Returns the incident matching edge from given node.
     2997    Edge matching(const Node& node) const {
     2998      return (*_matching)[node];
     2999    }
     3000
     3001    /// \brief Returns the mate of the node.
     3002    ///
     3003    /// Returns the adjancent node in a mathcing edge.
     3004    Node mate(const Node& node) const {
     3005      return _ugraph.target((*_matching)[node]);
     3006    }
     3007
     3008    /// @}
     3009
     3010    /// \name Dual solution
     3011    /// Functions for get the dual solution.
     3012   
     3013    /// @{
     3014
     3015    /// \brief Returns the value of the dual solution.
     3016    ///
     3017    /// Returns the value of the dual solution. It should be equal to
     3018    /// the primal value scaled by \ref dualScale "dual scale".
     3019    Value dualValue() const {
     3020      Value sum = 0;
     3021      for (NodeIt n(_ugraph); n != INVALID; ++n) {
     3022        sum += nodeValue(n);
     3023      }
     3024      for (int i = 0; i < blossomNum(); ++i) {
     3025        sum += blossomValue(i) * (blossomSize(i) / 2);
     3026      }
     3027      return sum;
     3028    }
     3029
     3030    /// \brief Returns the value of the node.
     3031    ///
     3032    /// Returns the the value of the node.
     3033    Value nodeValue(const Node& n) const {
     3034      return (*_node_potential)[n];
     3035    }
     3036
     3037    /// \brief Returns the number of the blossoms in the basis.
     3038    ///
     3039    /// Returns the number of the blossoms in the basis.
     3040    /// \see BlossomIt
     3041    int blossomNum() const {
     3042      return _blossom_potential.size();
     3043    }
     3044
     3045
     3046    /// \brief Returns the number of the nodes in the blossom.
     3047    ///
     3048    /// Returns the number of the nodes in the blossom.
     3049    int blossomSize(int k) const {
     3050      return _blossom_potential[k].end - _blossom_potential[k].begin;
     3051    }
     3052
     3053    /// \brief Returns the value of the blossom.
     3054    ///
     3055    /// Returns the the value of the blossom.
     3056    /// \see BlossomIt
     3057    Value blossomValue(int k) const {
     3058      return _blossom_potential[k].value;
     3059    }
     3060
     3061    /// \brief Lemon iterator for get the items of the blossom.
     3062    ///
     3063    /// Lemon iterator for get the nodes of the blossom. This class
     3064    /// provides a common style lemon iterator which gives back a
     3065    /// subset of the nodes.
     3066    class BlossomIt {
     3067    public:
     3068
     3069      /// \brief Constructor.
     3070      ///
     3071      /// Constructor for get the nodes of the variable.
     3072      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
     3073        : _algorithm(&algorithm)
     3074      {
     3075        _index = _algorithm->_blossom_potential[variable].begin;
     3076        _last = _algorithm->_blossom_potential[variable].end;
     3077      }
     3078
     3079      /// \brief Invalid constructor.
     3080      ///
     3081      /// Invalid constructor.
     3082      BlossomIt(Invalid) : _index(-1) {}
     3083
     3084      /// \brief Conversion to node.
     3085      ///
     3086      /// Conversion to node.
     3087      operator Node() const {
     3088        return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
     3089      }
     3090
     3091      /// \brief Increment operator.
     3092      ///
     3093      /// Increment operator.
     3094      BlossomIt& operator++() {
     3095        ++_index;
     3096        if (_index == _last) {
     3097          _index = -1;
     3098        }
     3099        return *this;
     3100      }
     3101
     3102      bool operator==(const BlossomIt& it) const {
     3103        return _index == it._index;
     3104      }
     3105      bool operator!=(const BlossomIt& it) const {
     3106        return _index != it._index;
     3107      }
     3108
     3109    private:
     3110      const MaxWeightedPerfectMatching* _algorithm;
     3111      int _last;
     3112      int _index;
     3113    };
     3114
     3115    /// @}
     3116   
     3117  };
     3118
    6213119 
    6223120} //END OF NAMESPACE LEMON
  • lemon/unionfind.h

    r2542 r2548  
    925925  };
    926926
     927  /// \ingroup auxdat
     928  ///
     929  /// \brief A \e Union-Find data structure implementation which
     930  /// is able to store a priority for each item and retrieve the minimum of
     931  /// each class.
     932  ///
     933  /// A \e Union-Find data structure implementation which is able to
     934  /// store a priority for each item and retrieve the minimum of each
     935  /// class. In addition, it supports the joining and splitting the
     936  /// components. If you don't need this feature then you makes
     937  /// better to use the \ref UnionFind class which is more efficient.
     938  ///
     939  /// The union-find data strcuture based on a (2, 16)-tree with a
     940  /// tournament minimum selection on the internal nodes. The insert
     941  /// operation takes O(1), the find, set, decrease and increase takes
     942  /// O(log(n)), where n is the number of nodes in the current
     943  /// component.  The complexity of join and split is O(log(n)*k),
     944  /// where n is the sum of the number of the nodes and k is the
     945  /// number of joined components or the number of the components
     946  /// after the split.
     947  ///
     948  /// \pre You need to add all the elements by the \ref insert()
     949  /// method.
     950  ///
     951  template <typename _Value, typename _ItemIntMap,
     952            typename _Comp = std::less<_Value> >
     953  class HeapUnionFind {
     954  public:
     955   
     956    typedef _Value Value;
     957    typedef typename _ItemIntMap::Key Item;
     958
     959    typedef _ItemIntMap ItemIntMap;
     960
     961    typedef _Comp Comp;
     962
     963  private:
     964
     965    static const int cmax = 3;
     966
     967    ItemIntMap& index;
     968
     969    struct ClassNode {
     970      int parent;
     971      int depth;
     972
     973      int left, right;
     974      int next, prev;
     975    };
     976
     977    int first_class;
     978    int first_free_class;
     979    std::vector<ClassNode> classes;
     980
     981    int newClass() {
     982      if (first_free_class < 0) {
     983        int id = classes.size();
     984        classes.push_back(ClassNode());
     985        return id;
     986      } else {
     987        int id = first_free_class;
     988        first_free_class = classes[id].next;
     989        return id;
     990      }
     991    }
     992
     993    void deleteClass(int id) {
     994      classes[id].next = first_free_class;
     995      first_free_class = id;
     996    }
     997
     998    struct ItemNode {
     999      int parent;
     1000      Item item;
     1001      Value prio;
     1002      int next, prev;
     1003      int left, right;
     1004      int size;
     1005    };
     1006
     1007    int first_free_node;
     1008    std::vector<ItemNode> nodes;
     1009
     1010    int newNode() {
     1011      if (first_free_node < 0) {
     1012        int id = nodes.size();
     1013        nodes.push_back(ItemNode());
     1014        return id;
     1015      } else {
     1016        int id = first_free_node;
     1017        first_free_node = nodes[id].next;
     1018        return id;
     1019      }
     1020    }
     1021
     1022    void deleteNode(int id) {
     1023      nodes[id].next = first_free_node;
     1024      first_free_node = id;
     1025    }
     1026
     1027    Comp comp;
     1028
     1029    int findClass(int id) const {
     1030      int kd = id;
     1031      while (kd >= 0) {
     1032        kd = nodes[kd].parent;
     1033      }
     1034      return ~kd;
     1035    }
     1036
     1037    int leftNode(int id) const {
     1038      int kd = ~(classes[id].parent);
     1039      for (int i = 0; i < classes[id].depth; ++i) {
     1040        kd = nodes[kd].left;
     1041      }
     1042      return kd;
     1043    }
     1044
     1045    int nextNode(int id) const {
     1046      int depth = 0;
     1047      while (id >= 0 && nodes[id].next == -1) {
     1048        id = nodes[id].parent;
     1049        ++depth;
     1050      }
     1051      if (id < 0) {
     1052        return -1;
     1053      }
     1054      id = nodes[id].next;
     1055      while (depth--) {
     1056        id = nodes[id].left;     
     1057      }
     1058      return id;
     1059    }
     1060
     1061
     1062    void setPrio(int id) {
     1063      int jd = nodes[id].left;
     1064      nodes[id].prio = nodes[jd].prio;
     1065      nodes[id].item = nodes[jd].item;
     1066      jd = nodes[jd].next;
     1067      while (jd != -1) {
     1068        if (comp(nodes[jd].prio, nodes[id].prio)) {
     1069          nodes[id].prio = nodes[jd].prio;
     1070          nodes[id].item = nodes[jd].item;
     1071        }
     1072        jd = nodes[jd].next;
     1073      }
     1074    }
     1075
     1076    void push(int id, int jd) {
     1077      nodes[id].size = 1;
     1078      nodes[id].left = nodes[id].right = jd;
     1079      nodes[jd].next = nodes[jd].prev = -1;
     1080      nodes[jd].parent = id;
     1081    }
     1082
     1083    void pushAfter(int id, int jd) {
     1084      int kd = nodes[id].parent;
     1085      if (nodes[id].next != -1) {
     1086        nodes[nodes[id].next].prev = jd;
     1087        if (kd >= 0) {
     1088          nodes[kd].size += 1;
     1089        }
     1090      } else {
     1091        if (kd >= 0) {
     1092          nodes[kd].right = jd;
     1093          nodes[kd].size += 1;
     1094        }
     1095      }
     1096      nodes[jd].next = nodes[id].next;
     1097      nodes[jd].prev = id;
     1098      nodes[id].next = jd;
     1099      nodes[jd].parent = kd;
     1100    }
     1101
     1102    void pushRight(int id, int jd) {
     1103      nodes[id].size += 1;
     1104      nodes[jd].prev = nodes[id].right;
     1105      nodes[jd].next = -1;
     1106      nodes[nodes[id].right].next = jd;
     1107      nodes[id].right = jd;
     1108      nodes[jd].parent = id;
     1109    }
     1110
     1111    void popRight(int id) {
     1112      nodes[id].size -= 1;
     1113      int jd = nodes[id].right;
     1114      nodes[nodes[jd].prev].next = -1;
     1115      nodes[id].right = nodes[jd].prev;
     1116    }
     1117
     1118    void splice(int id, int jd) {
     1119      nodes[id].size += nodes[jd].size;
     1120      nodes[nodes[id].right].next = nodes[jd].left;
     1121      nodes[nodes[jd].left].prev = nodes[id].right;
     1122      int kd = nodes[jd].left;
     1123      while (kd != -1) {
     1124        nodes[kd].parent = id;
     1125        kd = nodes[kd].next;
     1126      }
     1127    }
     1128
     1129    void split(int id, int jd) {
     1130      int kd = nodes[id].parent;
     1131      nodes[kd].right = nodes[id].prev;
     1132      nodes[nodes[id].prev].next = -1;
     1133     
     1134      nodes[jd].left = id;
     1135      nodes[id].prev = -1;
     1136      int num = 0;
     1137      while (id != -1) {
     1138        nodes[id].parent = jd;
     1139        nodes[jd].right = id;
     1140        id = nodes[id].next;
     1141        ++num;
     1142      }     
     1143      nodes[kd].size -= num;
     1144      nodes[jd].size = num;
     1145    }
     1146
     1147    void pushLeft(int id, int jd) {
     1148      nodes[id].size += 1;
     1149      nodes[jd].next = nodes[id].left;
     1150      nodes[jd].prev = -1;
     1151      nodes[nodes[id].left].prev = jd;
     1152      nodes[id].left = jd;
     1153      nodes[jd].parent = id;
     1154    }
     1155
     1156    void popLeft(int id) {
     1157      nodes[id].size -= 1;
     1158      int jd = nodes[id].left;
     1159      nodes[nodes[jd].next].prev = -1;
     1160      nodes[id].left = nodes[jd].next;
     1161    }
     1162
     1163    void repairLeft(int id) {
     1164      int jd = ~(classes[id].parent);
     1165      while (nodes[jd].left != -1) {
     1166        int kd = nodes[jd].left;
     1167        if (nodes[jd].size == 1) {
     1168          if (nodes[jd].parent < 0) {
     1169            classes[id].parent = ~kd;
     1170            classes[id].depth -= 1;
     1171            nodes[kd].parent = ~id;
     1172            deleteNode(jd);
     1173            jd = kd;
     1174          } else {
     1175            int pd = nodes[jd].parent;
     1176            if (nodes[nodes[jd].next].size < cmax) {
     1177              pushLeft(nodes[jd].next, nodes[jd].left);
     1178              if (less(nodes[jd].left, nodes[jd].next)) {
     1179                nodes[nodes[jd].next].prio = nodes[nodes[jd].left].prio;
     1180                nodes[nodes[jd].next].item = nodes[nodes[jd].left].item;
     1181              }
     1182              popLeft(pd);
     1183              deleteNode(jd);
     1184              jd = pd;
     1185            } else {
     1186              int ld = nodes[nodes[jd].next].left;
     1187              popLeft(nodes[jd].next);
     1188              pushRight(jd, ld);
     1189              if (less(ld, nodes[jd].left)) {
     1190                nodes[jd].item = nodes[ld].item;
     1191                nodes[jd].prio = nodes[jd].prio;
     1192              }
     1193              if (nodes[nodes[jd].next].item == nodes[ld].item) {
     1194                setPrio(nodes[jd].next);
     1195              }
     1196              jd = nodes[jd].left;
     1197            }
     1198          }
     1199        } else {
     1200          jd = nodes[jd].left;
     1201        }
     1202      }
     1203    }   
     1204
     1205    void repairRight(int id) {
     1206      int jd = ~(classes[id].parent);
     1207      while (nodes[jd].right != -1) {
     1208        int kd = nodes[jd].right;
     1209        if (nodes[jd].size == 1) {
     1210          if (nodes[jd].parent < 0) {
     1211            classes[id].parent = ~kd;
     1212            classes[id].depth -= 1;
     1213            nodes[kd].parent = ~id;
     1214            deleteNode(jd);
     1215            jd = kd;
     1216          } else {
     1217            int pd = nodes[jd].parent;
     1218            if (nodes[nodes[jd].prev].size < cmax) {
     1219              pushRight(nodes[jd].prev, nodes[jd].right);
     1220              if (less(nodes[jd].right, nodes[jd].prev)) {
     1221                nodes[nodes[jd].prev].prio = nodes[nodes[jd].right].prio;
     1222                nodes[nodes[jd].prev].item = nodes[nodes[jd].right].item;
     1223              }
     1224              popRight(pd);
     1225              deleteNode(jd);
     1226              jd = pd;
     1227            } else {
     1228              int ld = nodes[nodes[jd].prev].right;
     1229              popRight(nodes[jd].prev);
     1230              pushLeft(jd, ld);
     1231              if (less(ld, nodes[jd].right)) {
     1232                nodes[jd].item = nodes[ld].item;
     1233                nodes[jd].prio = nodes[jd].prio;
     1234              }
     1235              if (nodes[nodes[jd].prev].item == nodes[ld].item) {
     1236                setPrio(nodes[jd].prev);
     1237              }
     1238              jd = nodes[jd].right;
     1239            }
     1240          }
     1241        } else {
     1242          jd = nodes[jd].right;
     1243        }
     1244      }
     1245    }
     1246
     1247
     1248    bool less(int id, int jd) const {
     1249      return comp(nodes[id].prio, nodes[jd].prio);
     1250    }
     1251
     1252    bool equal(int id, int jd) const {
     1253      return !less(id, jd) && !less(jd, id);
     1254    }
     1255
     1256
     1257  public:
     1258
     1259    /// \brief Returns true when the given class is alive.
     1260    ///
     1261    /// Returns true when the given class is alive, ie. the class is
     1262    /// not nested into other class.
     1263    bool alive(int cls) const {
     1264      return classes[cls].parent < 0;
     1265    }
     1266
     1267    /// \brief Returns true when the given class is trivial.
     1268    ///
     1269    /// Returns true when the given class is trivial, ie. the class
     1270    /// contains just one item directly.
     1271    bool trivial(int cls) const {
     1272      return classes[cls].left == -1;
     1273    }
     1274
     1275    /// \brief Constructs the union-find.
     1276    ///
     1277    /// Constructs the union-find. 
     1278    /// \brief _index The index map of the union-find. The data
     1279    /// structure uses internally for store references.
     1280    HeapUnionFind(ItemIntMap& _index)
     1281      : index(_index), first_class(-1),
     1282        first_free_class(-1), first_free_node(-1) {}
     1283
     1284    /// \brief Insert a new node into a new component.
     1285    ///
     1286    /// Insert a new node into a new component.
     1287    /// \param item The item of the new node.
     1288    /// \param prio The priority of the new node.
     1289    /// \return The class id of the one-item-heap.
     1290    int insert(const Item& item, const Value& prio) {
     1291      int id = newNode();
     1292      nodes[id].item = item;
     1293      nodes[id].prio = prio;
     1294      nodes[id].size = 0;
     1295
     1296      nodes[id].prev = -1;
     1297      nodes[id].next = -1;
     1298
     1299      nodes[id].left = -1;
     1300      nodes[id].right = -1;
     1301
     1302      nodes[id].item = item;
     1303      index[item] = id;
     1304     
     1305      int class_id = newClass();
     1306      classes[class_id].parent = ~id;
     1307      classes[class_id].depth = 0;
     1308
     1309      classes[class_id].left = -1;
     1310      classes[class_id].right = -1;
     1311     
     1312      if (first_class != -1) {
     1313        classes[first_class].prev = class_id;
     1314      }
     1315      classes[class_id].next = first_class;
     1316      classes[class_id].prev = -1;
     1317      first_class = class_id;
     1318
     1319      nodes[id].parent = ~class_id;
     1320     
     1321      return class_id;
     1322    }
     1323
     1324    /// \brief The class of the item.
     1325    ///
     1326    /// \return The alive class id of the item, which is not nested into
     1327    /// other classes.
     1328    ///
     1329    /// The time complexity is O(log(n)).
     1330    int find(const Item& item) const {
     1331      return findClass(index[item]);
     1332    }
     1333   
     1334    /// \brief Joins the classes.
     1335    ///
     1336    /// The current function joins the given classes. The parameter is
     1337    /// an STL range which should be contains valid class ids. The
     1338    /// time complexity is O(log(n)*k) where n is the overall number
     1339    /// of the joined nodes and k is the number of classes.
     1340    /// \return The class of the joined classes.
     1341    /// \pre The range should contain at least two class ids.
     1342    template <typename Iterator>
     1343    int join(Iterator begin, Iterator end) {
     1344      std::vector<int> cs;
     1345      for (Iterator it = begin; it != end; ++it) {
     1346        cs.push_back(*it);
     1347      }
     1348
     1349      int class_id = newClass();
     1350      { // creation union-find
     1351
     1352        if (first_class != -1) {
     1353          classes[first_class].prev = class_id;
     1354        }
     1355        classes[class_id].next = first_class;
     1356        classes[class_id].prev = -1;
     1357        first_class = class_id;
     1358
     1359        classes[class_id].depth = classes[cs[0]].depth;
     1360        classes[class_id].parent = classes[cs[0]].parent;
     1361        nodes[~(classes[class_id].parent)].parent = ~class_id;
     1362
     1363        int l = cs[0];
     1364
     1365        classes[class_id].left = l;
     1366        classes[class_id].right = l;
     1367
     1368        if (classes[l].next != -1) {
     1369          classes[classes[l].next].prev = classes[l].prev;
     1370        }
     1371        classes[classes[l].prev].next = classes[l].next;
     1372       
     1373        classes[l].prev = -1;
     1374        classes[l].next = -1;
     1375
     1376        classes[l].depth = leftNode(l);
     1377        classes[l].parent = class_id;
     1378       
     1379      }
     1380
     1381      { // merging of heap
     1382        int l = class_id;
     1383        for (int ci = 1; ci < int(cs.size()); ++ci) {
     1384          int r = cs[ci];
     1385          int rln = leftNode(r);
     1386          if (classes[l].depth > classes[r].depth) {
     1387            int id = ~(classes[l].parent);
     1388            for (int i = classes[r].depth + 1; i < classes[l].depth; ++i) {
     1389              id = nodes[id].right;
     1390            }
     1391            while (id >= 0 && nodes[id].size == cmax) {
     1392              int new_id = newNode();
     1393              int right_id = nodes[id].right;
     1394
     1395              popRight(id);
     1396              if (nodes[id].item == nodes[right_id].item) {
     1397                setPrio(id);
     1398              }
     1399              push(new_id, right_id);
     1400              pushRight(new_id, ~(classes[r].parent));
     1401              setPrio(new_id);
     1402
     1403              id = nodes[id].parent;
     1404              classes[r].parent = ~new_id;
     1405            }
     1406            if (id < 0) {
     1407              int new_parent = newNode();
     1408              nodes[new_parent].next = -1;
     1409              nodes[new_parent].prev = -1;
     1410              nodes[new_parent].parent = ~l;
     1411
     1412              push(new_parent, ~(classes[l].parent));
     1413              pushRight(new_parent, ~(classes[r].parent));
     1414              setPrio(new_parent);
     1415
     1416              classes[l].parent = ~new_parent;
     1417              classes[l].depth += 1;
     1418            } else {
     1419              pushRight(id, ~(classes[r].parent));
     1420              while (id >= 0 && less(~(classes[r].parent), id)) {
     1421                nodes[id].prio = nodes[~(classes[r].parent)].prio;
     1422                nodes[id].item = nodes[~(classes[r].parent)].item;
     1423                id = nodes[id].parent;
     1424              }
     1425            }
     1426          } else if (classes[r].depth > classes[l].depth) {
     1427            int id = ~(classes[r].parent);
     1428            for (int i = classes[l].depth + 1; i < classes[r].depth; ++i) {
     1429              id = nodes[id].left;
     1430            }
     1431            while (id >= 0 && nodes[id].size == cmax) {
     1432              int new_id = newNode();
     1433              int left_id = nodes[id].left;
     1434
     1435              popLeft(id);
     1436              if (nodes[id].prio == nodes[left_id].prio) {
     1437                setPrio(id);
     1438              }
     1439              push(new_id, left_id);
     1440              pushLeft(new_id, ~(classes[l].parent));
     1441              setPrio(new_id);
     1442
     1443              id = nodes[id].parent;
     1444              classes[l].parent = ~new_id;
     1445
     1446            }
     1447            if (id < 0) {
     1448              int new_parent = newNode();
     1449              nodes[new_parent].next = -1;
     1450              nodes[new_parent].prev = -1;
     1451              nodes[new_parent].parent = ~l;
     1452
     1453              push(new_parent, ~(classes[r].parent));
     1454              pushLeft(new_parent, ~(classes[l].parent));
     1455              setPrio(new_parent);
     1456             
     1457              classes[r].parent = ~new_parent;
     1458              classes[r].depth += 1;
     1459            } else {
     1460              pushLeft(id, ~(classes[l].parent));
     1461              while (id >= 0 && less(~(classes[l].parent), id)) {
     1462                nodes[id].prio = nodes[~(classes[l].parent)].prio;
     1463                nodes[id].item = nodes[~(classes[l].parent)].item;
     1464                id = nodes[id].parent;
     1465              }
     1466            }
     1467            nodes[~(classes[r].parent)].parent = ~l;
     1468            classes[l].parent = classes[r].parent;
     1469            classes[l].depth = classes[r].depth;
     1470          } else {
     1471            if (classes[l].depth != 0 &&
     1472                nodes[~(classes[l].parent)].size +
     1473                nodes[~(classes[r].parent)].size <= cmax) {
     1474              splice(~(classes[l].parent), ~(classes[r].parent));
     1475              deleteNode(~(classes[r].parent));
     1476              if (less(~(classes[r].parent), ~(classes[l].parent))) {
     1477                nodes[~(classes[l].parent)].prio =
     1478                  nodes[~(classes[r].parent)].prio;
     1479                nodes[~(classes[l].parent)].item =
     1480                  nodes[~(classes[r].parent)].item;
     1481              }
     1482            } else {
     1483              int new_parent = newNode();
     1484              nodes[new_parent].next = nodes[new_parent].prev = -1;
     1485              push(new_parent, ~(classes[l].parent));
     1486              pushRight(new_parent, ~(classes[r].parent));
     1487              setPrio(new_parent);
     1488           
     1489              classes[l].parent = ~new_parent;
     1490              classes[l].depth += 1;
     1491              nodes[new_parent].parent = ~l;
     1492            }
     1493          }
     1494          if (classes[r].next != -1) {
     1495            classes[classes[r].next].prev = classes[r].prev;
     1496          }
     1497          classes[classes[r].prev].next = classes[r].next;
     1498
     1499          classes[r].prev = classes[l].right;
     1500          classes[classes[l].right].next = r;
     1501          classes[l].right = r;
     1502          classes[r].parent = l;
     1503
     1504          classes[r].next = -1;
     1505          classes[r].depth = rln;
     1506        }
     1507      }
     1508      return class_id;
     1509    }
     1510
     1511    /// \brief Split the class to subclasses.
     1512    ///
     1513    /// The current function splits the given class. The join, which
     1514    /// made the current class, stored a reference to the
     1515    /// subclasses. The \c splitClass() member restores the classes
     1516    /// and creates the heaps. The parameter is an STL output iterator
     1517    /// which will be filled with the subclass ids. The time
     1518    /// complexity is O(log(n)*k) where n is the overall number of
     1519    /// nodes in the splitted classes and k is the number of the
     1520    /// classes.
     1521    template <typename Iterator>
     1522    void split(int cls, Iterator out) {
     1523      std::vector<int> cs;
     1524      { // splitting union-find
     1525        int id = cls;
     1526        int l = classes[id].left;
     1527
     1528        classes[l].parent = classes[id].parent;
     1529        classes[l].depth = classes[id].depth;
     1530
     1531        nodes[~(classes[l].parent)].parent = ~l;
     1532
     1533        *out++ = l;
     1534
     1535        while (l != -1) {
     1536          cs.push_back(l);
     1537          l = classes[l].next;
     1538        }
     1539
     1540        classes[classes[id].right].next = first_class;
     1541        classes[first_class].prev = classes[id].right;
     1542        first_class = classes[id].left;
     1543       
     1544        if (classes[id].next != -1) {
     1545          classes[classes[id].next].prev = classes[id].prev;
     1546        }
     1547        classes[classes[id].prev].next = classes[id].next;
     1548       
     1549        deleteClass(id);
     1550      }
     1551
     1552      {
     1553        for (int i = 1; i < int(cs.size()); ++i) {
     1554          int l = classes[cs[i]].depth;
     1555          while (nodes[nodes[l].parent].left == l) {
     1556            l = nodes[l].parent;
     1557          }
     1558          int r = l;   
     1559          while (nodes[l].parent >= 0) {
     1560            l = nodes[l].parent;
     1561            int new_node = newNode();
     1562
     1563            nodes[new_node].prev = -1;
     1564            nodes[new_node].next = -1;
     1565
     1566            split(r, new_node);
     1567            pushAfter(l, new_node);
     1568            setPrio(l);
     1569            setPrio(new_node);
     1570            r = new_node;
     1571          }
     1572          classes[cs[i]].parent = ~r;
     1573          classes[cs[i]].depth = classes[~(nodes[l].parent)].depth;
     1574          nodes[r].parent = ~cs[i];
     1575
     1576          nodes[l].next = -1;
     1577          nodes[r].prev = -1;
     1578
     1579          repairRight(~(nodes[l].parent));
     1580          repairLeft(cs[i]);
     1581         
     1582          *out++ = cs[i];
     1583        }
     1584      }
     1585    }
     1586
     1587    /// \brief Gives back the priority of the current item.
     1588    ///
     1589    /// \return Gives back the priority of the current item.
     1590    const Value& operator[](const Item& item) const {
     1591      return nodes[index[item]].prio;
     1592    }
     1593
     1594    /// \brief Sets the priority of the current item.
     1595    ///
     1596    /// Sets the priority of the current item.
     1597    void set(const Item& item, const Value& prio) {
     1598      if (comp(prio, nodes[index[item]].prio)) {
     1599        decrease(item, prio);
     1600      } else if (!comp(prio, nodes[index[item]].prio)) {
     1601        increase(item, prio);
     1602      }
     1603    }
     1604     
     1605    /// \brief Increase the priority of the current item.
     1606    ///
     1607    /// Increase the priority of the current item.
     1608    void increase(const Item& item, const Value& prio) {
     1609      int id = index[item];
     1610      int kd = nodes[id].parent;
     1611      nodes[id].prio = prio;
     1612      while (kd >= 0 && nodes[kd].item == item) {
     1613        setPrio(kd);
     1614        kd = nodes[kd].parent;
     1615      }
     1616    }
     1617
     1618    /// \brief Increase the priority of the current item.
     1619    ///
     1620    /// Increase the priority of the current item.
     1621    void decrease(const Item& item, const Value& prio) {
     1622      int id = index[item];
     1623      int kd = nodes[id].parent;
     1624      nodes[id].prio = prio;
     1625      while (kd >= 0 && less(id, kd)) {
     1626        nodes[kd].prio = prio;
     1627        nodes[kd].item = item;
     1628        kd = nodes[kd].parent;
     1629      }
     1630    }
     1631   
     1632    /// \brief Gives back the minimum priority of the class.
     1633    ///
     1634    /// \return Gives back the minimum priority of the class.
     1635    const Value& classPrio(int cls) const {
     1636      return nodes[~(classes[cls].parent)].prio;
     1637    }
     1638
     1639    /// \brief Gives back the minimum priority item of the class.
     1640    ///
     1641    /// \return Gives back the minimum priority item of the class.
     1642    const Item& classTop(int cls) const {
     1643      return nodes[~(classes[cls].parent)].item;
     1644    }
     1645
     1646    /// \brief Gives back a representant item of the class.
     1647    ///
     1648    /// The representant is indpendent from the priorities of the
     1649    /// items.
     1650    /// \return Gives back a representant item of the class.
     1651    const Item& classRep(int id) const {
     1652      int parent = classes[id].parent;
     1653      return nodes[parent >= 0 ? classes[id].depth : leftNode(id)].item;
     1654    }
     1655
     1656    /// \brief Lemon style iterator for the items of a class.
     1657    ///
     1658    /// ClassIt is a lemon style iterator for the components. It iterates
     1659    /// on the items of a class. By example if you want to iterate on
     1660    /// each items of each classes then you may write the next code.
     1661    ///\code
     1662    /// for (ClassIt cit(huf); cit != INVALID; ++cit) {
     1663    ///   std::cout << "Class: ";
     1664    ///   for (ItemIt iit(huf, cit); iit != INVALID; ++iit) {
     1665    ///     std::cout << toString(iit) << ' ' << std::endl;
     1666    ///   }
     1667    ///   std::cout << std::endl;
     1668    /// }
     1669    ///\endcode
     1670    class ItemIt {
     1671    private:
     1672
     1673      const HeapUnionFind* _huf;
     1674      int _id, _lid;
     1675     
     1676    public:
     1677
     1678      /// \brief Default constructor
     1679      ///
     1680      /// Default constructor
     1681      ItemIt() {}
     1682
     1683      ItemIt(const HeapUnionFind& huf, int cls) : _huf(&huf) {
     1684        int id = cls;
     1685        int parent = _huf->classes[id].parent;
     1686        if (parent >= 0) {
     1687          _id = _huf->classes[id].depth;
     1688          if (_huf->classes[id].next != -1) {
     1689            _lid = _huf->classes[_huf->classes[id].next].depth;
     1690          } else {
     1691            _lid = -1;
     1692          }
     1693        } else {
     1694          _id = _huf->leftNode(id);
     1695          _lid = -1;
     1696        }
     1697      }
     1698     
     1699      /// \brief Increment operator
     1700      ///
     1701      /// It steps to the next item in the class.
     1702      ItemIt& operator++() {
     1703        _id = _huf->nextNode(_id);
     1704        return *this;
     1705      }
     1706
     1707      /// \brief Conversion operator
     1708      ///
     1709      /// It converts the iterator to the current item.
     1710      operator const Item&() const {
     1711        return _huf->nodes[_id].item;
     1712      }
     1713     
     1714      /// \brief Equality operator
     1715      ///
     1716      /// Equality operator
     1717      bool operator==(const ItemIt& i) {
     1718        return i._id == _id;
     1719      }
     1720
     1721      /// \brief Inequality operator
     1722      ///
     1723      /// Inequality operator
     1724      bool operator!=(const ItemIt& i) {
     1725        return i._id != _id;
     1726      }
     1727
     1728      /// \brief Equality operator
     1729      ///
     1730      /// Equality operator
     1731      bool operator==(Invalid) {
     1732        return _id == _lid;
     1733      }
     1734
     1735      /// \brief Inequality operator
     1736      ///
     1737      /// Inequality operator
     1738      bool operator!=(Invalid) {
     1739        return _id != _lid;
     1740      }
     1741     
     1742    };
     1743
     1744    /// \brief Class iterator
     1745    ///
     1746    /// The iterator stores
     1747    class ClassIt {
     1748    private:
     1749
     1750      const HeapUnionFind* _huf;
     1751      int _id;
     1752
     1753    public:
     1754
     1755      ClassIt(const HeapUnionFind& huf)
     1756        : _huf(&huf), _id(huf.first_class) {}
     1757
     1758      ClassIt(const HeapUnionFind& huf, int cls)
     1759        : _huf(&huf), _id(huf.classes[cls].left) {}
     1760
     1761      ClassIt(Invalid) : _huf(0), _id(-1) {}
     1762     
     1763      const ClassIt& operator++() {
     1764        _id = _huf->classes[_id].next;
     1765        return *this;
     1766      }
     1767
     1768      /// \brief Equality operator
     1769      ///
     1770      /// Equality operator
     1771      bool operator==(const ClassIt& i) {
     1772        return i._id == _id;
     1773      }
     1774
     1775      /// \brief Inequality operator
     1776      ///
     1777      /// Inequality operator
     1778      bool operator!=(const ClassIt& i) {
     1779        return i._id != _id;
     1780      }     
     1781     
     1782      operator int() const {
     1783        return _id;
     1784      }
     1785           
     1786    };
     1787
     1788  };
     1789
    9271790  //! @}
    9281791
Note: See TracChangeset for help on using the changeset viewer.