COIN-OR::LEMON - Graph Library

Changeset 2548:a3ba22ebccc6 in lemon-0.x for lemon


Ignore:
Timestamp:
12/28/07 12:00:51 (12 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?

Location:
lemon
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • 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.