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

Ignore:
Timestamp:
12/28/07 12:00:51 (12 years ago)
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
 r2505 #define LEMON_MAX_MATCHING_H #include #include #include #include #include #include #include ///\ingroup matching ///\ingroup matching /// ///\brief Edmonds' alternating forest maximum matching algorithm. /// }; /// \ingroup matching /// /// \brief Weighted matching in general undirected graphs /// /// This class provides an efficient implementation of Edmond's /// maximum weighted matching algorithm. The implementation is based /// on extensive use of priority queues and provides /// \f$O(nm\log(n))\f$ time complexity. /// /// The maximum weighted matching problem is to find undirected /// edges in the graph with maximum overall weight and no two of /// them shares their endpoints. The problem can be formulated with /// the next linear program: /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f] ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f] /// \f[x_e \ge 0\quad \forall e\in E\f] /// \f[\max \sum_{e\in E}x_ew_e\f] /// where \f$\delta(X)\f$ is the set of edges incident to a node in /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of /// the nodes. /// /// The algorithm calculates an optimal matching and a proof of the /// optimality. The solution of the dual problem can be used to check /// the result of the algorithm. The dual linear problem is the next: /// \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] /// \f[y_u \ge 0 \quad \forall u \in V\f] /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f] /// /// The algorithm can be executed with \c run() or the \c init() and /// then the \c start() member functions. After it the matching can /// be asked with \c matching() or mate() functions. The dual /// solution can be get with \c nodeValue(), \c blossomNum() and \c /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt /// "BlossomIt" nested class which is able to iterate on the nodes /// of a blossom. If the value type is integral then the dual /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4". /// /// \author Balazs Dezso template > class MaxWeightedMatching { public: typedef _UGraph UGraph; typedef _WeightMap WeightMap; typedef typename WeightMap::Value Value; /// \brief Scaling factor for dual solution /// /// Scaling factor for dual solution, it is equal to 4 or 1 /// according to the value type. static const int dualScale = std::numeric_limits::is_integer ? 4 : 1; typedef typename UGraph::template NodeMap MatchingMap; private: UGRAPH_TYPEDEFS(typename UGraph); typedef typename UGraph::template NodeMap NodePotential; typedef std::vector BlossomNodeList; struct BlossomVariable { int begin, end; Value value; BlossomVariable(int _begin, int _end, Value _value) : begin(_begin), end(_end), value(_value) {} }; typedef std::vector BlossomPotential; const UGraph& _ugraph; const WeightMap& _weight; MatchingMap* _matching; NodePotential* _node_potential; BlossomPotential _blossom_potential; BlossomNodeList _blossom_node_list; int _node_num; int _blossom_num; typedef typename UGraph::template NodeMap NodeIntMap; typedef typename UGraph::template EdgeMap EdgeIntMap; typedef typename UGraph::template UEdgeMap UEdgeIntMap; typedef IntegerMap IntIntMap; enum Status { EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2 }; typedef HeapUnionFind BlossomSet; struct BlossomData { int tree; Status status; Edge pred, next; Value pot, offset; Node base; }; NodeIntMap *_blossom_index; BlossomSet *_blossom_set; IntegerMap* _blossom_data; NodeIntMap *_node_index; EdgeIntMap *_node_heap_index; struct NodeData { NodeData(EdgeIntMap& node_heap_index) : heap(node_heap_index) {} int blossom; Value pot; BinHeap heap; std::map heap_index; int tree; }; IntegerMap* _node_data; typedef ExtendFindEnum TreeSet; IntIntMap *_tree_set_index; TreeSet *_tree_set; NodeIntMap *_delta1_index; BinHeap *_delta1; IntIntMap *_delta2_index; BinHeap *_delta2; UEdgeIntMap *_delta3_index; BinHeap *_delta3; IntIntMap *_delta4_index; BinHeap *_delta4; Value _delta_sum; void createStructures() { _node_num = countNodes(_ugraph); _blossom_num = _node_num * 3 / 2; if (!_matching) { _matching = new MatchingMap(_ugraph); } if (!_node_potential) { _node_potential = new NodePotential(_ugraph); } if (!_blossom_set) { _blossom_index = new NodeIntMap(_ugraph); _blossom_set = new BlossomSet(*_blossom_index); _blossom_data = new IntegerMap(_blossom_num); } if (!_node_index) { _node_index = new NodeIntMap(_ugraph); _node_heap_index = new EdgeIntMap(_ugraph); _node_data = new IntegerMap(_node_num, NodeData(*_node_heap_index)); } if (!_tree_set) { _tree_set_index = new IntIntMap(_blossom_num); _tree_set = new TreeSet(*_tree_set_index); } if (!_delta1) { _delta1_index = new NodeIntMap(_ugraph); _delta1 = new BinHeap(*_delta1_index); } if (!_delta2) { _delta2_index = new IntIntMap(_blossom_num); _delta2 = new BinHeap(*_delta2_index); } if (!_delta3) { _delta3_index = new UEdgeIntMap(_ugraph); _delta3 = new BinHeap(*_delta3_index); } if (!_delta4) { _delta4_index = new IntIntMap(_blossom_num); _delta4 = new BinHeap(*_delta4_index); } } void destroyStructures() { _node_num = countNodes(_ugraph); _blossom_num = _node_num * 3 / 2; if (_matching) { delete _matching; } if (_node_potential) { delete _node_potential; } if (_blossom_set) { delete _blossom_index; delete _blossom_set; delete _blossom_data; } if (_node_index) { delete _node_index; delete _node_heap_index; delete _node_data; } if (_tree_set) { delete _tree_set_index; delete _tree_set; } if (_delta1) { delete _delta1_index; delete _delta1; } if (_delta2) { delete _delta2_index; delete _delta2; } if (_delta3) { delete _delta3_index; delete _delta3; } if (_delta4) { delete _delta4_index; delete _delta4; } } void matchedToEven(int blossom, int tree) { if (_delta2->state(blossom) == _delta2->IN_HEAP) { _delta2->erase(blossom); } if (!_blossom_set->trivial(blossom)) { (*_blossom_data)[blossom].pot -= 2 * (_delta_sum - (*_blossom_data)[blossom].offset); } for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); n != INVALID; ++n) { _blossom_set->increase(n, std::numeric_limits::max()); int ni = (*_node_index)[n]; (*_node_data)[ni].heap.clear(); (*_node_data)[ni].heap_index.clear(); (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset; _delta1->push(n, (*_node_data)[ni].pot); for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { Node v = _ugraph.source(e); int vb = _blossom_set->find(v); int vi = (*_node_index)[v]; Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - dualScale * _weight[e]; if ((*_blossom_data)[vb].status == EVEN) { if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { _delta3->push(e, rw / 2); } } else if ((*_blossom_data)[vb].status == UNMATCHED) { if (_delta3->state(e) != _delta3->IN_HEAP) { _delta3->push(e, rw); } } else { typename std::map::iterator it = (*_node_data)[vi].heap_index.find(tree); if (it != (*_node_data)[vi].heap_index.end()) { if ((*_node_data)[vi].heap[it->second] > rw) { (*_node_data)[vi].heap.replace(it->second, e); (*_node_data)[vi].heap.decrease(e, rw); it->second = e; } } else { (*_node_data)[vi].heap.push(e, rw); (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); } if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); if ((*_blossom_data)[vb].status == MATCHED) { if (_delta2->state(vb) != _delta2->IN_HEAP) { _delta2->push(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset){ _delta2->decrease(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); } } } } } } (*_blossom_data)[blossom].offset = 0; } void matchedToOdd(int blossom) { if (_delta2->state(blossom) == _delta2->IN_HEAP) { _delta2->erase(blossom); } (*_blossom_data)[blossom].offset += _delta_sum; if (!_blossom_set->trivial(blossom)) { _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + (*_blossom_data)[blossom].offset); } } void evenToMatched(int blossom, int tree) { if (!_blossom_set->trivial(blossom)) { (*_blossom_data)[blossom].pot += 2 * _delta_sum; } for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); n != INVALID; ++n) { int ni = (*_node_index)[n]; (*_node_data)[ni].pot -= _delta_sum; _delta1->erase(n); for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { Node v = _ugraph.source(e); int vb = _blossom_set->find(v); int vi = (*_node_index)[v]; Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - dualScale * _weight[e]; if (vb == blossom) { if (_delta3->state(e) == _delta3->IN_HEAP) { _delta3->erase(e); } } else if ((*_blossom_data)[vb].status == EVEN) { if (_delta3->state(e) == _delta3->IN_HEAP) { _delta3->erase(e); } int vt = _tree_set->find(vb); if (vt != tree) { Edge r = _ugraph.oppositeEdge(e); typename std::map::iterator it = (*_node_data)[ni].heap_index.find(vt); if (it != (*_node_data)[ni].heap_index.end()) { if ((*_node_data)[ni].heap[it->second] > rw) { (*_node_data)[ni].heap.replace(it->second, r); (*_node_data)[ni].heap.decrease(r, rw); it->second = r; } } else { (*_node_data)[ni].heap.push(r, rw); (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); } if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); if (_delta2->state(blossom) != _delta2->IN_HEAP) { _delta2->push(blossom, _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset); } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset){ _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset); } } } } else if ((*_blossom_data)[vb].status == UNMATCHED) { if (_delta3->state(e) == _delta3->IN_HEAP) { _delta3->erase(e); } } else { typename std::map::iterator it = (*_node_data)[vi].heap_index.find(tree); if (it != (*_node_data)[vi].heap_index.end()) { (*_node_data)[vi].heap.erase(it->second); (*_node_data)[vi].heap_index.erase(it); if ((*_node_data)[vi].heap.empty()) { _blossom_set->increase(v, std::numeric_limits::max()); } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); } if ((*_blossom_data)[vb].status == MATCHED) { if (_blossom_set->classPrio(vb) == std::numeric_limits::max()) { _delta2->erase(vb); } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset) { _delta2->increase(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); } } } } } } } void oddToMatched(int blossom) { (*_blossom_data)[blossom].offset -= _delta_sum; if (_blossom_set->classPrio(blossom) != std::numeric_limits::max()) { _delta2->push(blossom, _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset); } if (!_blossom_set->trivial(blossom)) { _delta4->erase(blossom); } } void oddToEven(int blossom, int tree) { if (!_blossom_set->trivial(blossom)) { _delta4->erase(blossom); (*_blossom_data)[blossom].pot -= 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); } for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); n != INVALID; ++n) { int ni = (*_node_index)[n]; _blossom_set->increase(n, std::numeric_limits::max()); (*_node_data)[ni].heap.clear(); (*_node_data)[ni].heap_index.clear(); (*_node_data)[ni].pot += 2 * _delta_sum - (*_blossom_data)[blossom].offset; _delta1->push(n, (*_node_data)[ni].pot); for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { Node v = _ugraph.source(e); int vb = _blossom_set->find(v); int vi = (*_node_index)[v]; Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - dualScale * _weight[e]; if ((*_blossom_data)[vb].status == EVEN) { if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { _delta3->push(e, rw / 2); } } else if ((*_blossom_data)[vb].status == UNMATCHED) { if (_delta3->state(e) != _delta3->IN_HEAP) { _delta3->push(e, rw); } } else { typename std::map::iterator it = (*_node_data)[vi].heap_index.find(tree); if (it != (*_node_data)[vi].heap_index.end()) { if ((*_node_data)[vi].heap[it->second] > rw) { (*_node_data)[vi].heap.replace(it->second, e); (*_node_data)[vi].heap.decrease(e, rw); it->second = e; } } else { (*_node_data)[vi].heap.push(e, rw); (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); } if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); if ((*_blossom_data)[vb].status == MATCHED) { if (_delta2->state(vb) != _delta2->IN_HEAP) { _delta2->push(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset) { _delta2->decrease(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); } } } } } } (*_blossom_data)[blossom].offset = 0; } void matchedToUnmatched(int blossom) { if (_delta2->state(blossom) == _delta2->IN_HEAP) { _delta2->erase(blossom); } for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); n != INVALID; ++n) { int ni = (*_node_index)[n]; _blossom_set->increase(n, std::numeric_limits::max()); (*_node_data)[ni].heap.clear(); (*_node_data)[ni].heap_index.clear(); for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) { Node v = _ugraph.target(e); int vb = _blossom_set->find(v); int vi = (*_node_index)[v]; Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - dualScale * _weight[e]; if ((*_blossom_data)[vb].status == EVEN) { if (_delta3->state(e) != _delta3->IN_HEAP) { _delta3->push(e, rw); } } } } } void unmatchedToMatched(int blossom) { for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); n != INVALID; ++n) { int ni = (*_node_index)[n]; for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { Node v = _ugraph.source(e); int vb = _blossom_set->find(v); int vi = (*_node_index)[v]; Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - dualScale * _weight[e]; if (vb == blossom) { if (_delta3->state(e) == _delta3->IN_HEAP) { _delta3->erase(e); } } else if ((*_blossom_data)[vb].status == EVEN) { if (_delta3->state(e) == _delta3->IN_HEAP) { _delta3->erase(e); } int vt = _tree_set->find(vb); Edge r = _ugraph.oppositeEdge(e); typename std::map::iterator it = (*_node_data)[ni].heap_index.find(vt); if (it != (*_node_data)[ni].heap_index.end()) { if ((*_node_data)[ni].heap[it->second] > rw) { (*_node_data)[ni].heap.replace(it->second, r); (*_node_data)[ni].heap.decrease(r, rw); it->second = r; } } else { (*_node_data)[ni].heap.push(r, rw); (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); } if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); if (_delta2->state(blossom) != _delta2->IN_HEAP) { _delta2->push(blossom, _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset); } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)- (*_blossom_data)[blossom].offset){ _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset); } } } else if ((*_blossom_data)[vb].status == UNMATCHED) { if (_delta3->state(e) == _delta3->IN_HEAP) { _delta3->erase(e); } } } } } void alternatePath(int even, int tree) { int odd; evenToMatched(even, tree); (*_blossom_data)[even].status = MATCHED; while ((*_blossom_data)[even].pred != INVALID) { odd = _blossom_set->find(_ugraph.target((*_blossom_data)[even].pred)); (*_blossom_data)[odd].status = MATCHED; oddToMatched(odd); (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred)); (*_blossom_data)[even].status = MATCHED; evenToMatched(even, tree); (*_blossom_data)[even].next = _ugraph.oppositeEdge((*_blossom_data)[odd].pred); } } void destroyTree(int tree) { for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { if ((*_blossom_data)[b].status == EVEN) { (*_blossom_data)[b].status = MATCHED; evenToMatched(b, tree); } else if ((*_blossom_data)[b].status == ODD) { (*_blossom_data)[b].status = MATCHED; oddToMatched(b); } } _tree_set->eraseClass(tree); } void unmatchNode(const Node& node) { int blossom = _blossom_set->find(node); int tree = _tree_set->find(blossom); alternatePath(blossom, tree); destroyTree(tree); (*_blossom_data)[blossom].status = UNMATCHED; (*_blossom_data)[blossom].base = node; matchedToUnmatched(blossom); } void augmentOnEdge(const UEdge& edge) { int left = _blossom_set->find(_ugraph.source(edge)); int right = _blossom_set->find(_ugraph.target(edge)); if ((*_blossom_data)[left].status == EVEN) { int left_tree = _tree_set->find(left); alternatePath(left, left_tree); destroyTree(left_tree); } else { (*_blossom_data)[left].status = MATCHED; unmatchedToMatched(left); } if ((*_blossom_data)[right].status == EVEN) { int right_tree = _tree_set->find(right); alternatePath(right, right_tree); destroyTree(right_tree); } else { (*_blossom_data)[right].status = MATCHED; unmatchedToMatched(right); } (*_blossom_data)[left].next = _ugraph.direct(edge, true); (*_blossom_data)[right].next = _ugraph.direct(edge, false); } void extendOnEdge(const Edge& edge) { int base = _blossom_set->find(_ugraph.target(edge)); int tree = _tree_set->find(base); int odd = _blossom_set->find(_ugraph.source(edge)); _tree_set->insert(odd, tree); (*_blossom_data)[odd].status = ODD; matchedToOdd(odd); (*_blossom_data)[odd].pred = edge; int even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].next)); (*_blossom_data)[even].pred = (*_blossom_data)[even].next; _tree_set->insert(even, tree); (*_blossom_data)[even].status = EVEN; matchedToEven(even, tree); } void shrinkOnEdge(const UEdge& uedge, int tree) { int nca = -1; std::vector left_path, right_path; { std::set left_set, right_set; int left = _blossom_set->find(_ugraph.source(uedge)); left_path.push_back(left); left_set.insert(left); int right = _blossom_set->find(_ugraph.target(uedge)); right_path.push_back(right); right_set.insert(right); while (true) { if ((*_blossom_data)[left].pred == INVALID) break; left = _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred)); left_path.push_back(left); left = _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred)); left_path.push_back(left); left_set.insert(left); if (right_set.find(left) != right_set.end()) { nca = left; break; } if ((*_blossom_data)[right].pred == INVALID) break; right = _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred)); right_path.push_back(right); right = _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred)); right_path.push_back(right); right_set.insert(right); if (left_set.find(right) != left_set.end()) { nca = right; break; } } if (nca == -1) { if ((*_blossom_data)[left].pred == INVALID) { nca = right; while (left_set.find(nca) == left_set.end()) { nca = _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred)); right_path.push_back(nca); nca = _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred)); right_path.push_back(nca); } } else { nca = left; while (right_set.find(nca) == right_set.end()) { nca = _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred)); left_path.push_back(nca); nca = _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred)); left_path.push_back(nca); } } } } std::vector subblossoms; Edge prev; prev = _ugraph.direct(uedge, true); for (int i = 0; left_path[i] != nca; i += 2) { subblossoms.push_back(left_path[i]); (*_blossom_data)[left_path[i]].next = prev; _tree_set->erase(left_path[i]); subblossoms.push_back(left_path[i + 1]); (*_blossom_data)[left_path[i + 1]].status = EVEN; oddToEven(left_path[i + 1], tree); _tree_set->erase(left_path[i + 1]); prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred); } int k = 0; while (right_path[k] != nca) ++k; subblossoms.push_back(nca); (*_blossom_data)[nca].next = prev; for (int i = k - 2; i >= 0; i -= 2) { subblossoms.push_back(right_path[i + 1]); (*_blossom_data)[right_path[i + 1]].status = EVEN; oddToEven(right_path[i + 1], tree); _tree_set->erase(right_path[i + 1]); (*_blossom_data)[right_path[i + 1]].next = (*_blossom_data)[right_path[i + 1]].pred; subblossoms.push_back(right_path[i]); _tree_set->erase(right_path[i]); } int surface = _blossom_set->join(subblossoms.begin(), subblossoms.end()); for (int i = 0; i < int(subblossoms.size()); ++i) { if (!_blossom_set->trivial(subblossoms[i])) { (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; } (*_blossom_data)[subblossoms[i]].status = MATCHED; } (*_blossom_data)[surface].pot = -2 * _delta_sum; (*_blossom_data)[surface].offset = 0; (*_blossom_data)[surface].status = EVEN; (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; _tree_set->insert(surface, tree); _tree_set->erase(nca); } void splitBlossom(int blossom) { Edge next = (*_blossom_data)[blossom].next; Edge pred = (*_blossom_data)[blossom].pred; int tree = _tree_set->find(blossom); (*_blossom_data)[blossom].status = MATCHED; oddToMatched(blossom); if (_delta2->state(blossom) == _delta2->IN_HEAP) { _delta2->erase(blossom); } std::vector subblossoms; _blossom_set->split(blossom, std::back_inserter(subblossoms)); Value offset = (*_blossom_data)[blossom].offset; int b = _blossom_set->find(_ugraph.source(pred)); int d = _blossom_set->find(_ugraph.source(next)); int ib, id; for (int i = 0; i < int(subblossoms.size()); ++i) { if (subblossoms[i] == b) ib = i; if (subblossoms[i] == d) id = i; (*_blossom_data)[subblossoms[i]].offset = offset; if (!_blossom_set->trivial(subblossoms[i])) { (*_blossom_data)[subblossoms[i]].pot -= 2 * offset; } if (_blossom_set->classPrio(subblossoms[i]) != std::numeric_limits::max()) { _delta2->push(subblossoms[i], _blossom_set->classPrio(subblossoms[i]) - (*_blossom_data)[subblossoms[i]].offset); } } if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) { for (int i = (id + 1) % subblossoms.size(); i != ib; i = (i + 2) % subblossoms.size()) { int sb = subblossoms[i]; int tb = subblossoms[(i + 1) % subblossoms.size()]; (*_blossom_data)[sb].next = _ugraph.oppositeEdge((*_blossom_data)[tb].next); } for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { int sb = subblossoms[i]; int tb = subblossoms[(i + 1) % subblossoms.size()]; int ub = subblossoms[(i + 2) % subblossoms.size()]; (*_blossom_data)[sb].status = ODD; matchedToOdd(sb); _tree_set->insert(sb, tree); (*_blossom_data)[sb].pred = pred; (*_blossom_data)[sb].next = _ugraph.oppositeEdge((*_blossom_data)[tb].next); pred = (*_blossom_data)[ub].next; (*_blossom_data)[tb].status = EVEN; matchedToEven(tb, tree); _tree_set->insert(tb, tree); (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; } (*_blossom_data)[subblossoms[id]].status = ODD; matchedToOdd(subblossoms[id]); _tree_set->insert(subblossoms[id], tree); (*_blossom_data)[subblossoms[id]].next = next; (*_blossom_data)[subblossoms[id]].pred = pred; } else { for (int i = (ib + 1) % subblossoms.size(); i != id; i = (i + 2) % subblossoms.size()) { int sb = subblossoms[i]; int tb = subblossoms[(i + 1) % subblossoms.size()]; (*_blossom_data)[sb].next = _ugraph.oppositeEdge((*_blossom_data)[tb].next); } for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { int sb = subblossoms[i]; int tb = subblossoms[(i + 1) % subblossoms.size()]; int ub = subblossoms[(i + 2) % subblossoms.size()]; (*_blossom_data)[sb].status = ODD; matchedToOdd(sb); _tree_set->insert(sb, tree); (*_blossom_data)[sb].next = next; (*_blossom_data)[sb].pred = _ugraph.oppositeEdge((*_blossom_data)[tb].next); (*_blossom_data)[tb].status = EVEN; matchedToEven(tb, tree); _tree_set->insert(tb, tree); (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next = _ugraph.oppositeEdge((*_blossom_data)[ub].next); next = (*_blossom_data)[ub].next; } (*_blossom_data)[subblossoms[ib]].status = ODD; matchedToOdd(subblossoms[ib]); _tree_set->insert(subblossoms[ib], tree); (*_blossom_data)[subblossoms[ib]].next = next; (*_blossom_data)[subblossoms[ib]].pred = pred; } _tree_set->erase(blossom); } void extractBlossom(int blossom, const Node& base, const Edge& matching) { if (_blossom_set->trivial(blossom)) { int bi = (*_node_index)[base]; Value pot = (*_node_data)[bi].pot; _matching->set(base, matching); _blossom_node_list.push_back(base); _node_potential->set(base, pot); } else { Value pot = (*_blossom_data)[blossom].pot; int bn = _blossom_node_list.size(); std::vector subblossoms; _blossom_set->split(blossom, std::back_inserter(subblossoms)); int b = _blossom_set->find(base); int ib = -1; for (int i = 0; i < int(subblossoms.size()); ++i) { if (subblossoms[i] == b) { ib = i; break; } } for (int i = 1; i < int(subblossoms.size()); i += 2) { int sb = subblossoms[(ib + i) % subblossoms.size()]; int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; Edge m = (*_blossom_data)[tb].next; extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m)); extractBlossom(tb, _ugraph.source(m), m); } extractBlossom(subblossoms[ib], base, matching); int en = _blossom_node_list.size(); _blossom_potential.push_back(BlossomVariable(bn, en, pot)); } } void extractMatching() { std::vector blossoms; for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { blossoms.push_back(c); } for (int i = 0; i < int(blossoms.size()); ++i) { if ((*_blossom_data)[blossoms[i]].status == MATCHED) { Value offset = (*_blossom_data)[blossoms[i]].offset; (*_blossom_data)[blossoms[i]].pot += 2 * offset; for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); n != INVALID; ++n) { (*_node_data)[(*_node_index)[n]].pot -= offset; } Edge matching = (*_blossom_data)[blossoms[i]].next; Node base = _ugraph.source(matching); extractBlossom(blossoms[i], base, matching); } else { Node base = (*_blossom_data)[blossoms[i]].base; extractBlossom(blossoms[i], base, INVALID); } } } public: /// \brief Constructor /// /// Constructor. MaxWeightedMatching(const UGraph& ugraph, const WeightMap& weight) : _ugraph(ugraph), _weight(weight), _matching(0), _node_potential(0), _blossom_potential(), _blossom_node_list(), _node_num(0), _blossom_num(0), _blossom_index(0), _blossom_set(0), _blossom_data(0), _node_index(0), _node_heap_index(0), _node_data(0), _tree_set_index(0), _tree_set(0), _delta1_index(0), _delta1(0), _delta2_index(0), _delta2(0), _delta3_index(0), _delta3(0), _delta4_index(0), _delta4(0), _delta_sum() {} ~MaxWeightedMatching() { destroyStructures(); } /// \name Execution control /// The simplest way to execute the algorithm is to use the member /// \c run() member function. ///@{ /// \brief Initialize the algorithm /// /// Initialize the algorithm void init() { createStructures(); for (EdgeIt e(_ugraph); e != INVALID; ++e) { _node_heap_index->set(e, BinHeap::PRE_HEAP); } for (NodeIt n(_ugraph); n != INVALID; ++n) { _delta1_index->set(n, _delta1->PRE_HEAP); } for (UEdgeIt e(_ugraph); e != INVALID; ++e) { _delta3_index->set(e, _delta3->PRE_HEAP); } for (int i = 0; i < _blossom_num; ++i) { _delta2_index->set(i, _delta2->PRE_HEAP); _delta4_index->set(i, _delta4->PRE_HEAP); } int index = 0; for (NodeIt n(_ugraph); n != INVALID; ++n) { Value max = 0; for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) { if (_ugraph.target(e) == n) continue; if ((dualScale * _weight[e]) / 2 > max) { max = (dualScale * _weight[e]) / 2; } } _node_index->set(n, index); (*_node_data)[index].pot = max; _delta1->push(n, max); int blossom = _blossom_set->insert(n, std::numeric_limits::max()); _tree_set->insert(blossom); (*_blossom_data)[blossom].status = EVEN; (*_blossom_data)[blossom].pred = INVALID; (*_blossom_data)[blossom].next = INVALID; (*_blossom_data)[blossom].pot = 0; (*_blossom_data)[blossom].offset = 0; ++index; } for (UEdgeIt e(_ugraph); e != INVALID; ++e) { int si = (*_node_index)[_ugraph.source(e)]; int ti = (*_node_index)[_ugraph.target(e)]; if (_ugraph.source(e) != _ugraph.target(e)) { _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - dualScale * _weight[e]) / 2); } } } /// \brief Starts the algorithm /// /// Starts the algorithm void start() { enum OpType { D1, D2, D3, D4 }; int unmatched = _node_num; while (unmatched > 0) { Value d1 = !_delta1->empty() ? _delta1->prio() : std::numeric_limits::max(); Value d2 = !_delta2->empty() ? _delta2->prio() : std::numeric_limits::max(); Value d3 = !_delta3->empty() ? _delta3->prio() : std::numeric_limits::max(); Value d4 = !_delta4->empty() ? _delta4->prio() : std::numeric_limits::max(); _delta_sum = d1; OpType ot = D1; if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } switch (ot) { case D1: { Node n = _delta1->top(); unmatchNode(n); --unmatched; } break; case D2: { int blossom = _delta2->top(); Node n = _blossom_set->classTop(blossom); Edge e = (*_node_data)[(*_node_index)[n]].heap.top(); extendOnEdge(e); } break; case D3: { UEdge e = _delta3->top(); int left_blossom = _blossom_set->find(_ugraph.source(e)); int right_blossom = _blossom_set->find(_ugraph.target(e)); if (left_blossom == right_blossom) { _delta3->pop(); } else { int left_tree; if ((*_blossom_data)[left_blossom].status == EVEN) { left_tree = _tree_set->find(left_blossom); } else { left_tree = -1; ++unmatched; } int right_tree; if ((*_blossom_data)[right_blossom].status == EVEN) { right_tree = _tree_set->find(right_blossom); } else { right_tree = -1; ++unmatched; } if (left_tree == right_tree) { shrinkOnEdge(e, left_tree); } else { augmentOnEdge(e); unmatched -= 2; } } } break; case D4: splitBlossom(_delta4->top()); break; } } extractMatching(); } /// \brief Runs %MaxWeightedMatching algorithm. /// /// This method runs the %MaxWeightedMatching algorithm. /// /// \note mwm.run() is just a shortcut of the following code. /// \code ///   mwm.init(); ///   mwm.start(); /// \endcode void run() { init(); start(); } /// @} /// \name Primal solution /// Functions for get the primal solution, ie. the matching. /// @{ /// \brief Returns the matching value. /// /// Returns the matching value. Value matchingValue() const { Value sum = 0; for (NodeIt n(_ugraph); n != INVALID; ++n) { if ((*_matching)[n] != INVALID) { sum += _weight[(*_matching)[n]]; } } return sum /= 2; } /// \brief Returns true when the edge is in the matching. /// /// Returns true when the edge is in the matching. bool matching(const UEdge& edge) const { return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true); } /// \brief Returns the incident matching edge. /// /// Returns the incident matching edge from given node. If the /// node is not matched then it gives back \c INVALID. Edge matching(const Node& node) const { return (*_matching)[node]; } /// \brief Returns the mate of the node. /// /// Returns the adjancent node in a mathcing edge. If the node is /// not matched then it gives back \c INVALID. Node mate(const Node& node) const { return (*_matching)[node] != INVALID ? _ugraph.target((*_matching)[node]) : INVALID; } /// @} /// \name Dual solution /// Functions for get the dual solution. /// @{ /// \brief Returns the value of the dual solution. /// /// Returns the value of the dual solution. It should be equal to /// the primal value scaled by \ref dualScale "dual scale". Value dualValue() const { Value sum = 0; for (NodeIt n(_ugraph); n != INVALID; ++n) { sum += nodeValue(n); } for (int i = 0; i < blossomNum(); ++i) { sum += blossomValue(i) * (blossomSize(i) / 2); } return sum; } /// \brief Returns the value of the node. /// /// Returns the the value of the node. Value nodeValue(const Node& n) const { return (*_node_potential)[n]; } /// \brief Returns the number of the blossoms in the basis. /// /// Returns the number of the blossoms in the basis. /// \see BlossomIt int blossomNum() const { return _blossom_potential.size(); } /// \brief Returns the number of the nodes in the blossom. /// /// Returns the number of the nodes in the blossom. int blossomSize(int k) const { return _blossom_potential[k].end - _blossom_potential[k].begin; } /// \brief Returns the value of the blossom. /// /// Returns the the value of the blossom. /// \see BlossomIt Value blossomValue(int k) const { return _blossom_potential[k].value; } /// \brief Lemon iterator for get the items of the blossom. /// /// Lemon iterator for get the nodes of the blossom. This class /// provides a common style lemon iterator which gives back a /// subset of the nodes. class BlossomIt { public: /// \brief Constructor. /// /// Constructor for get the nodes of the variable. BlossomIt(const MaxWeightedMatching& algorithm, int variable) : _algorithm(&algorithm) { _index = _algorithm->_blossom_potential[variable].begin; _last = _algorithm->_blossom_potential[variable].end; } /// \brief Invalid constructor. /// /// Invalid constructor. BlossomIt(Invalid) : _index(-1) {} /// \brief Conversion to node. /// /// Conversion to node. operator Node() const { return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID; } /// \brief Increment operator. /// /// Increment operator. BlossomIt& operator++() { ++_index; if (_index == _last) { _index = -1; } return *this; } bool operator==(const BlossomIt& it) const { return _index == it._index; } bool operator!=(const BlossomIt& it) const { return _index != it._index; } private: const MaxWeightedMatching* _algorithm; int _last; int _index; }; /// @} }; /// \ingroup matching /// /// \brief Weighted perfect matching in general undirected graphs /// /// This class provides an efficient implementation of Edmond's /// maximum weighted perfecr matching algorithm. The implementation /// is based on extensive use of priority queues and provides /// \f$O(nm\log(n))\f$ time complexity. /// /// The maximum weighted matching problem is to find undirected /// edges in the graph with maximum overall weight and no two of /// them shares their endpoints and covers all nodes. The problem /// can be formulated with the next linear program: /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f] ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f] /// \f[x_e \ge 0\quad \forall e\in E\f] /// \f[\max \sum_{e\in E}x_ew_e\f] /// where \f$\delta(X)\f$ is the set of edges incident to a node in /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of /// the nodes. /// /// The algorithm calculates an optimal matching and a proof of the /// optimality. The solution of the dual problem can be used to check /// the result of the algorithm. The dual linear problem is the next: /// \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] /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f] /// /// The algorithm can be executed with \c run() or the \c init() and /// then the \c start() member functions. After it the matching can /// be asked with \c matching() or mate() functions. The dual /// solution can be get with \c nodeValue(), \c blossomNum() and \c /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt /// "BlossomIt" nested class which is able to iterate on the nodes /// of a blossom. If the value type is integral then the dual /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4". /// /// \author Balazs Dezso template > class MaxWeightedPerfectMatching { public: typedef _UGraph UGraph; typedef _WeightMap WeightMap; typedef typename WeightMap::Value Value; /// \brief Scaling factor for dual solution /// /// Scaling factor for dual solution, it is equal to 4 or 1 /// according to the value type. static const int dualScale = std::numeric_limits::is_integer ? 4 : 1; typedef typename UGraph::template NodeMap MatchingMap; private: UGRAPH_TYPEDEFS(typename UGraph); typedef typename UGraph::template NodeMap NodePotential; typedef std::vector BlossomNodeList; struct BlossomVariable { int begin, end; Value value; BlossomVariable(int _begin, int _end, Value _value) : begin(_begin), end(_end), value(_value) {} }; typedef std::vector BlossomPotential; const UGraph& _ugraph; const WeightMap& _weight; MatchingMap* _matching; NodePotential* _node_potential; BlossomPotential _blossom_potential; BlossomNodeList _blossom_node_list; int _node_num; int _blossom_num; typedef typename UGraph::template NodeMap NodeIntMap; typedef typename UGraph::template EdgeMap EdgeIntMap; typedef typename UGraph::template UEdgeMap UEdgeIntMap; typedef IntegerMap IntIntMap; enum Status { EVEN = -1, MATCHED = 0, ODD = 1 }; typedef HeapUnionFind BlossomSet; struct BlossomData { int tree; Status status; Edge pred, next; Value pot, offset; }; NodeIntMap *_blossom_index; BlossomSet *_blossom_set; IntegerMap* _blossom_data; NodeIntMap *_node_index; EdgeIntMap *_node_heap_index; struct NodeData { NodeData(EdgeIntMap& node_heap_index) : heap(node_heap_index) {} int blossom; Value pot; BinHeap heap; std::map heap_index; int tree; }; IntegerMap* _node_data; typedef ExtendFindEnum TreeSet; IntIntMap *_tree_set_index; TreeSet *_tree_set; IntIntMap *_delta2_index; BinHeap *_delta2; UEdgeIntMap *_delta3_index; BinHeap *_delta3; IntIntMap *_delta4_index; BinHeap *_delta4; Value _delta_sum; void createStructures() { _node_num = countNodes(_ugraph); _blossom_num = _node_num * 3 / 2; if (!_matching) { _matching = new MatchingMap(_ugraph); } if (!_node_potential) { _node_potential = new NodePotential(_ugraph); } if (!_blossom_set) { _blossom_index = new NodeIntMap(_ugraph); _blossom_set = new BlossomSet(*_blossom_index); _blossom_data = new IntegerMap(_blossom_num); } if (!_node_index) { _node_index = new NodeIntMap(_ugraph); _node_heap_index = new EdgeIntMap(_ugraph); _node_data = new IntegerMap(_node_num, NodeData(*_node_heap_index)); } if (!_tree_set) { _tree_set_index = new IntIntMap(_blossom_num); _tree_set = new TreeSet(*_tree_set_index); } if (!_delta2) { _delta2_index = new IntIntMap(_blossom_num); _delta2 = new BinHeap(*_delta2_index); } if (!_delta3) { _delta3_index = new UEdgeIntMap(_ugraph); _delta3 = new BinHeap(*_delta3_index); } if (!_delta4) { _delta4_index = new IntIntMap(_blossom_num); _delta4 = new BinHeap(*_delta4_index); } } void destroyStructures() { _node_num = countNodes(_ugraph); _blossom_num = _node_num * 3 / 2; if (_matching) { delete _matching; } if (_node_potential) { delete _node_potential; } if (_blossom_set) { delete _blossom_index; delete _blossom_set; delete _blossom_data; } if (_node_index) { delete _node_index; delete _node_heap_index; delete _node_data; } if (_tree_set) { delete _tree_set_index; delete _tree_set; } if (_delta2) { delete _delta2_index; delete _delta2; } if (_delta3) { delete _delta3_index; delete _delta3; } if (_delta4) { delete _delta4_index; delete _delta4; } } void matchedToEven(int blossom, int tree) { if (_delta2->state(blossom) == _delta2->IN_HEAP) { _delta2->erase(blossom); } if (!_blossom_set->trivial(blossom)) { (*_blossom_data)[blossom].pot -= 2 * (_delta_sum - (*_blossom_data)[blossom].offset); } for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); n != INVALID; ++n) { _blossom_set->increase(n, std::numeric_limits::max()); int ni = (*_node_index)[n]; (*_node_data)[ni].heap.clear(); (*_node_data)[ni].heap_index.clear(); (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset; for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { Node v = _ugraph.source(e); int vb = _blossom_set->find(v); int vi = (*_node_index)[v]; Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - dualScale * _weight[e]; if ((*_blossom_data)[vb].status == EVEN) { if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { _delta3->push(e, rw / 2); } } else { typename std::map::iterator it = (*_node_data)[vi].heap_index.find(tree); if (it != (*_node_data)[vi].heap_index.end()) { if ((*_node_data)[vi].heap[it->second] > rw) { (*_node_data)[vi].heap.replace(it->second, e); (*_node_data)[vi].heap.decrease(e, rw); it->second = e; } } else { (*_node_data)[vi].heap.push(e, rw); (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); } if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); if ((*_blossom_data)[vb].status == MATCHED) { if (_delta2->state(vb) != _delta2->IN_HEAP) { _delta2->push(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset){ _delta2->decrease(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); } } } } } } (*_blossom_data)[blossom].offset = 0; } void matchedToOdd(int blossom) { if (_delta2->state(blossom) == _delta2->IN_HEAP) { _delta2->erase(blossom); } (*_blossom_data)[blossom].offset += _delta_sum; if (!_blossom_set->trivial(blossom)) { _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + (*_blossom_data)[blossom].offset); } } void evenToMatched(int blossom, int tree) { if (!_blossom_set->trivial(blossom)) { (*_blossom_data)[blossom].pot += 2 * _delta_sum; } for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); n != INVALID; ++n) { int ni = (*_node_index)[n]; (*_node_data)[ni].pot -= _delta_sum; for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { Node v = _ugraph.source(e); int vb = _blossom_set->find(v); int vi = (*_node_index)[v]; Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - dualScale * _weight[e]; if (vb == blossom) { if (_delta3->state(e) == _delta3->IN_HEAP) { _delta3->erase(e); } } else if ((*_blossom_data)[vb].status == EVEN) { if (_delta3->state(e) == _delta3->IN_HEAP) { _delta3->erase(e); } int vt = _tree_set->find(vb); if (vt != tree) { Edge r = _ugraph.oppositeEdge(e); typename std::map::iterator it = (*_node_data)[ni].heap_index.find(vt); if (it != (*_node_data)[ni].heap_index.end()) { if ((*_node_data)[ni].heap[it->second] > rw) { (*_node_data)[ni].heap.replace(it->second, r); (*_node_data)[ni].heap.decrease(r, rw); it->second = r; } } else { (*_node_data)[ni].heap.push(r, rw); (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); } if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); if (_delta2->state(blossom) != _delta2->IN_HEAP) { _delta2->push(blossom, _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset); } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset){ _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset); } } } } else { typename std::map::iterator it = (*_node_data)[vi].heap_index.find(tree); if (it != (*_node_data)[vi].heap_index.end()) { (*_node_data)[vi].heap.erase(it->second); (*_node_data)[vi].heap_index.erase(it); if ((*_node_data)[vi].heap.empty()) { _blossom_set->increase(v, std::numeric_limits::max()); } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); } if ((*_blossom_data)[vb].status == MATCHED) { if (_blossom_set->classPrio(vb) == std::numeric_limits::max()) { _delta2->erase(vb); } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset) { _delta2->increase(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); } } } } } } } void oddToMatched(int blossom) { (*_blossom_data)[blossom].offset -= _delta_sum; if (_blossom_set->classPrio(blossom) != std::numeric_limits::max()) { _delta2->push(blossom, _blossom_set->classPrio(blossom) - (*_blossom_data)[blossom].offset); } if (!_blossom_set->trivial(blossom)) { _delta4->erase(blossom); } } void oddToEven(int blossom, int tree) { if (!_blossom_set->trivial(blossom)) { _delta4->erase(blossom); (*_blossom_data)[blossom].pot -= 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); } for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); n != INVALID; ++n) { int ni = (*_node_index)[n]; _blossom_set->increase(n, std::numeric_limits::max()); (*_node_data)[ni].heap.clear(); (*_node_data)[ni].heap_index.clear(); (*_node_data)[ni].pot += 2 * _delta_sum - (*_blossom_data)[blossom].offset; for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) { Node v = _ugraph.source(e); int vb = _blossom_set->find(v); int vi = (*_node_index)[v]; Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - dualScale * _weight[e]; if ((*_blossom_data)[vb].status == EVEN) { if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { _delta3->push(e, rw / 2); } } else { typename std::map::iterator it = (*_node_data)[vi].heap_index.find(tree); if (it != (*_node_data)[vi].heap_index.end()) { if ((*_node_data)[vi].heap[it->second] > rw) { (*_node_data)[vi].heap.replace(it->second, e); (*_node_data)[vi].heap.decrease(e, rw); it->second = e; } } else { (*_node_data)[vi].heap.push(e, rw); (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); } if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); if ((*_blossom_data)[vb].status == MATCHED) { if (_delta2->state(vb) != _delta2->IN_HEAP) { _delta2->push(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset) { _delta2->decrease(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); } } } } } } (*_blossom_data)[blossom].offset = 0; } void alternatePath(int even, int tree) { int odd; evenToMatched(even, tree); (*_blossom_data)[even].status = MATCHED; while ((*_blossom_data)[even].pred != INVALID) { odd = _blossom_set->find(_ugraph.target((*_blossom_data)[even].pred)); (*_blossom_data)[odd].status = MATCHED; oddToMatched(odd); (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred)); (*_blossom_data)[even].status = MATCHED; evenToMatched(even, tree); (*_blossom_data)[even].next = _ugraph.oppositeEdge((*_blossom_data)[odd].pred); } } void destroyTree(int tree) { for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { if ((*_blossom_data)[b].status == EVEN) { (*_blossom_data)[b].status = MATCHED; evenToMatched(b, tree); } else if ((*_blossom_data)[b].status == ODD) { (*_blossom_data)[b].status = MATCHED; oddToMatched(b); } } _tree_set->eraseClass(tree); } void augmentOnEdge(const UEdge& edge) { int left = _blossom_set->find(_ugraph.source(edge)); int right = _blossom_set->find(_ugraph.target(edge)); int left_tree = _tree_set->find(left); alternatePath(left, left_tree); destroyTree(left_tree); int right_tree = _tree_set->find(right); alternatePath(right, right_tree); destroyTree(right_tree); (*_blossom_data)[left].next = _ugraph.direct(edge, true); (*_blossom_data)[right].next = _ugraph.direct(edge, false); } void extendOnEdge(const Edge& edge) { int base = _blossom_set->find(_ugraph.target(edge)); int tree = _tree_set->find(base); int odd = _blossom_set->find(_ugraph.source(edge)); _tree_set->insert(odd, tree); (*_blossom_data)[odd].status = ODD; matchedToOdd(odd); (*_blossom_data)[odd].pred = edge; int even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].next)); (*_blossom_data)[even].pred = (*_blossom_data)[even].next; _tree_set->insert(even, tree); (*_blossom_data)[even].status = EVEN; matchedToEven(even, tree); } void shrinkOnEdge(const UEdge& uedge, int tree) { int nca = -1; std::vector left_path, right_path; { std::set left_set, right_set; int left = _blossom_set->find(_ugraph.source(uedge)); left_path.push_back(left); left_set.insert(left); int right = _blossom_set->find(_ugraph.target(uedge)); right_path.push_back(right); right_set.insert(right); while (true) { if ((*_blossom_data)[left].pred == INVALID) break; left = _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred)); left_path.push_back(left); left = _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred)); left_path.push_back(left); left_set.insert(left); if (right_set.find(left) != right_set.end()) { nca = left; break; } if ((*_blossom_data)[right].pred == INVALID) break; right = _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred)); right_path.push_back(right); right = _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred)); right_path.push_back(right); right_set.insert(right); if (left_set.find(right) != left_set.end()) { nca = right; break; } } if (nca == -1) { if ((*_blossom_data)[left].pred == INVALID) { nca = right; while (left_set.find(nca) == left_set.end()) { nca = _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred)); right_path.push_back(nca); nca = _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred)); right_path.push_back(nca); } } else { nca = left; while (right_set.find(nca) == right_set.end()) { nca = _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred)); left_path.push_back(nca); nca = _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred)); left_path.push_back(nca); } } } } std::vector subblossoms; Edge prev; prev = _ugraph.direct(uedge, true); for (int i = 0; left_path[i] != nca; i += 2) { subblossoms.push_back(left_path[i]); (*_blossom_data)[left_path[i]].next = prev; _tree_set->erase(left_path[i]); subblossoms.push_back(left_path[i + 1]); (*_blossom_data)[left_path[i + 1]].status = EVEN; oddToEven(left_path[i + 1], tree); _tree_set->erase(left_path[i + 1]); prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred); } int k = 0; while (right_path[k] != nca) ++k; subblossoms.push_back(nca); (*_blossom_data)[nca].next = prev; for (int i = k - 2; i >= 0; i -= 2) { subblossoms.push_back(right_path[i + 1]); (*_blossom_data)[right_path[i + 1]].status = EVEN; oddToEven(right_path[i + 1], tree); _tree_set->erase(right_path[i + 1]); (*_blossom_data)[right_path[i + 1]].next = (*_blossom_data)[right_path[i + 1]].pred; subblossoms.push_back(right_path[i]); _tree_set->erase(right_path[i]); } int surface = _blossom_set->join(subblossoms.begin(), subblossoms.end()); for (int i = 0; i < int(subblossoms.size()); ++i) { if (!_blossom_set->trivial(subblossoms[i])) { (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; } (*_blossom_data)[subblossoms[i]].status = MATCHED; } (*_blossom_data)[surface].pot = -2 * _delta_sum; (*_blossom_data)[surface].offset = 0; (*_blossom_data)[surface].status = EVEN; (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; _tree_set->insert(surface, tree); _tree_set->erase(nca); } void splitBlossom(int blossom) { Edge next = (*_blossom_data)[blossom].next; Edge pred = (*_blossom_data)[blossom].pred; int tree = _tree_set->find(blossom); (*_blossom_data)[blossom].status = MATCHED; oddToMatched(blossom); if (_delta2->state(blossom) == _delta2->IN_HEAP) { _delta2->erase(blossom); } std::vector subblossoms; _blossom_set->split(blossom, std::back_inserter(subblossoms)); Value offset = (*_blossom_data)[blossom].offset; int b = _blossom_set->find(_ugraph.source(pred)); int d = _blossom_set->find(_ugraph.source(next)); int ib, id; for (int i = 0; i < int(subblossoms.size()); ++i) { if (subblossoms[i] == b) ib = i; if (subblossoms[i] == d) id = i; (*_blossom_data)[subblossoms[i]].offset = offset; if (!_blossom_set->trivial(subblossoms[i])) { (*_blossom_data)[subblossoms[i]].pot -= 2 * offset; } if (_blossom_set->classPrio(subblossoms[i]) != std::numeric_limits::max()) { _delta2->push(subblossoms[i], _blossom_set->classPrio(subblossoms[i]) - (*_blossom_data)[subblossoms[i]].offset); } } if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) { for (int i = (id + 1) % subblossoms.size(); i != ib; i = (i + 2) % subblossoms.size()) { int sb = subblossoms[i]; int tb = subblossoms[(i + 1) % subblossoms.size()]; (*_blossom_data)[sb].next = _ugraph.oppositeEdge((*_blossom_data)[tb].next); } for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { int sb = subblossoms[i]; int tb = subblossoms[(i + 1) % subblossoms.size()]; int ub = subblossoms[(i + 2) % subblossoms.size()]; (*_blossom_data)[sb].status = ODD; matchedToOdd(sb); _tree_set->insert(sb, tree); (*_blossom_data)[sb].pred = pred; (*_blossom_data)[sb].next = _ugraph.oppositeEdge((*_blossom_data)[tb].next); pred = (*_blossom_data)[ub].next; (*_blossom_data)[tb].status = EVEN; matchedToEven(tb, tree); _tree_set->insert(tb, tree); (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; } (*_blossom_data)[subblossoms[id]].status = ODD; matchedToOdd(subblossoms[id]); _tree_set->insert(subblossoms[id], tree); (*_blossom_data)[subblossoms[id]].next = next; (*_blossom_data)[subblossoms[id]].pred = pred; } else { for (int i = (ib + 1) % subblossoms.size(); i != id; i = (i + 2) % subblossoms.size()) { int sb = subblossoms[i]; int tb = subblossoms[(i + 1) % subblossoms.size()]; (*_blossom_data)[sb].next = _ugraph.oppositeEdge((*_blossom_data)[tb].next); } for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { int sb = subblossoms[i]; int tb = subblossoms[(i + 1) % subblossoms.size()]; int ub = subblossoms[(i + 2) % subblossoms.size()]; (*_blossom_data)[sb].status = ODD; matchedToOdd(sb); _tree_set->insert(sb, tree); (*_blossom_data)[sb].next = next; (*_blossom_data)[sb].pred = _ugraph.oppositeEdge((*_blossom_data)[tb].next); (*_blossom_data)[tb].status = EVEN; matchedToEven(tb, tree); _tree_set->insert(tb, tree); (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next = _ugraph.oppositeEdge((*_blossom_data)[ub].next); next = (*_blossom_data)[ub].next; } (*_blossom_data)[subblossoms[ib]].status = ODD; matchedToOdd(subblossoms[ib]); _tree_set->insert(subblossoms[ib], tree); (*_blossom_data)[subblossoms[ib]].next = next; (*_blossom_data)[subblossoms[ib]].pred = pred; } _tree_set->erase(blossom); } void extractBlossom(int blossom, const Node& base, const Edge& matching) { if (_blossom_set->trivial(blossom)) { int bi = (*_node_index)[base]; Value pot = (*_node_data)[bi].pot; _matching->set(base, matching); _blossom_node_list.push_back(base); _node_potential->set(base, pot); } else { Value pot = (*_blossom_data)[blossom].pot; int bn = _blossom_node_list.size(); std::vector subblossoms; _blossom_set->split(blossom, std::back_inserter(subblossoms)); int b = _blossom_set->find(base); int ib = -1; for (int i = 0; i < int(subblossoms.size()); ++i) { if (subblossoms[i] == b) { ib = i; break; } } for (int i = 1; i < int(subblossoms.size()); i += 2) { int sb = subblossoms[(ib + i) % subblossoms.size()]; int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; Edge m = (*_blossom_data)[tb].next; extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m)); extractBlossom(tb, _ugraph.source(m), m); } extractBlossom(subblossoms[ib], base, matching); int en = _blossom_node_list.size(); _blossom_potential.push_back(BlossomVariable(bn, en, pot)); } } void extractMatching() { std::vector blossoms; for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { blossoms.push_back(c); } for (int i = 0; i < int(blossoms.size()); ++i) { Value offset = (*_blossom_data)[blossoms[i]].offset; (*_blossom_data)[blossoms[i]].pot += 2 * offset; for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); n != INVALID; ++n) { (*_node_data)[(*_node_index)[n]].pot -= offset; } Edge matching = (*_blossom_data)[blossoms[i]].next; Node base = _ugraph.source(matching); extractBlossom(blossoms[i], base, matching); } } public: /// \brief Constructor /// /// Constructor. MaxWeightedPerfectMatching(const UGraph& ugraph, const WeightMap& weight) : _ugraph(ugraph), _weight(weight), _matching(0), _node_potential(0), _blossom_potential(), _blossom_node_list(), _node_num(0), _blossom_num(0), _blossom_index(0), _blossom_set(0), _blossom_data(0), _node_index(0), _node_heap_index(0), _node_data(0), _tree_set_index(0), _tree_set(0), _delta2_index(0), _delta2(0), _delta3_index(0), _delta3(0), _delta4_index(0), _delta4(0), _delta_sum() {} ~MaxWeightedPerfectMatching() { destroyStructures(); } /// \name Execution control /// The simplest way to execute the algorithm is to use the member /// \c run() member function. ///@{ /// \brief Initialize the algorithm /// /// Initialize the algorithm void init() { createStructures(); for (EdgeIt e(_ugraph); e != INVALID; ++e) { _node_heap_index->set(e, BinHeap::PRE_HEAP); } for (UEdgeIt e(_ugraph); e != INVALID; ++e) { _delta3_index->set(e, _delta3->PRE_HEAP); } for (int i = 0; i < _blossom_num; ++i) { _delta2_index->set(i, _delta2->PRE_HEAP); _delta4_index->set(i, _delta4->PRE_HEAP); } int index = 0; for (NodeIt n(_ugraph); n != INVALID; ++n) { Value max = std::numeric_limits::min(); for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) { if (_ugraph.target(e) == n) continue; if ((dualScale * _weight[e]) / 2 > max) { max = (dualScale * _weight[e]) / 2; } } _node_index->set(n, index); (*_node_data)[index].pot = max; int blossom = _blossom_set->insert(n, std::numeric_limits::max()); _tree_set->insert(blossom); (*_blossom_data)[blossom].status = EVEN; (*_blossom_data)[blossom].pred = INVALID; (*_blossom_data)[blossom].next = INVALID; (*_blossom_data)[blossom].pot = 0; (*_blossom_data)[blossom].offset = 0; ++index; } for (UEdgeIt e(_ugraph); e != INVALID; ++e) { int si = (*_node_index)[_ugraph.source(e)]; int ti = (*_node_index)[_ugraph.target(e)]; if (_ugraph.source(e) != _ugraph.target(e)) { _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - dualScale * _weight[e]) / 2); } } } /// \brief Starts the algorithm /// /// Starts the algorithm bool start() { enum OpType { D2, D3, D4 }; int unmatched = _node_num; while (unmatched > 0) { Value d2 = !_delta2->empty() ? _delta2->prio() : std::numeric_limits::max(); Value d3 = !_delta3->empty() ? _delta3->prio() : std::numeric_limits::max(); Value d4 = !_delta4->empty() ? _delta4->prio() : std::numeric_limits::max(); _delta_sum = d2; OpType ot = D2; if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } if (_delta_sum == std::numeric_limits::max()) { return false; } switch (ot) { case D2: { int blossom = _delta2->top(); Node n = _blossom_set->classTop(blossom); Edge e = (*_node_data)[(*_node_index)[n]].heap.top(); extendOnEdge(e); } break; case D3: { UEdge e = _delta3->top(); int left_blossom = _blossom_set->find(_ugraph.source(e)); int right_blossom = _blossom_set->find(_ugraph.target(e)); if (left_blossom == right_blossom) { _delta3->pop(); } else { int left_tree = _tree_set->find(left_blossom); int right_tree = _tree_set->find(right_blossom); if (left_tree == right_tree) { shrinkOnEdge(e, left_tree); } else { augmentOnEdge(e); unmatched -= 2; } } } break; case D4: splitBlossom(_delta4->top()); break; } } extractMatching(); return true; } /// \brief Runs %MaxWeightedPerfectMatching algorithm. /// /// This method runs the %MaxWeightedPerfectMatching algorithm. /// /// \note mwm.run() is just a shortcut of the following code. /// \code ///   mwm.init(); ///   mwm.start(); /// \endcode bool run() { init(); return start(); } /// @} /// \name Primal solution /// Functions for get the primal solution, ie. the matching. /// @{ /// \brief Returns the matching value. /// /// Returns the matching value. Value matchingValue() const { Value sum = 0; for (NodeIt n(_ugraph); n != INVALID; ++n) { if ((*_matching)[n] != INVALID) { sum += _weight[(*_matching)[n]]; } } return sum /= 2; } /// \brief Returns true when the edge is in the matching. /// /// Returns true when the edge is in the matching. bool matching(const UEdge& edge) const { return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true); } /// \brief Returns the incident matching edge. /// /// Returns the incident matching edge from given node. Edge matching(const Node& node) const { return (*_matching)[node]; } /// \brief Returns the mate of the node. /// /// Returns the adjancent node in a mathcing edge. Node mate(const Node& node) const { return _ugraph.target((*_matching)[node]); } /// @} /// \name Dual solution /// Functions for get the dual solution. /// @{ /// \brief Returns the value of the dual solution. /// /// Returns the value of the dual solution. It should be equal to /// the primal value scaled by \ref dualScale "dual scale". Value dualValue() const { Value sum = 0; for (NodeIt n(_ugraph); n != INVALID; ++n) { sum += nodeValue(n); } for (int i = 0; i < blossomNum(); ++i) { sum += blossomValue(i) * (blossomSize(i) / 2); } return sum; } /// \brief Returns the value of the node. /// /// Returns the the value of the node. Value nodeValue(const Node& n) const { return (*_node_potential)[n]; } /// \brief Returns the number of the blossoms in the basis. /// /// Returns the number of the blossoms in the basis. /// \see BlossomIt int blossomNum() const { return _blossom_potential.size(); } /// \brief Returns the number of the nodes in the blossom. /// /// Returns the number of the nodes in the blossom. int blossomSize(int k) const { return _blossom_potential[k].end - _blossom_potential[k].begin; } /// \brief Returns the value of the blossom. /// /// Returns the the value of the blossom. /// \see BlossomIt Value blossomValue(int k) const { return _blossom_potential[k].value; } /// \brief Lemon iterator for get the items of the blossom. /// /// Lemon iterator for get the nodes of the blossom. This class /// provides a common style lemon iterator which gives back a /// subset of the nodes. class BlossomIt { public: /// \brief Constructor. /// /// Constructor for get the nodes of the variable. BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) : _algorithm(&algorithm) { _index = _algorithm->_blossom_potential[variable].begin; _last = _algorithm->_blossom_potential[variable].end; } /// \brief Invalid constructor. /// /// Invalid constructor. BlossomIt(Invalid) : _index(-1) {} /// \brief Conversion to node. /// /// Conversion to node. operator Node() const { return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID; } /// \brief Increment operator. /// /// Increment operator. BlossomIt& operator++() { ++_index; if (_index == _last) { _index = -1; } return *this; } bool operator==(const BlossomIt& it) const { return _index == it._index; } bool operator!=(const BlossomIt& it) const { return _index != it._index; } private: const MaxWeightedPerfectMatching* _algorithm; int _last; int _index; }; /// @} }; } //END OF NAMESPACE LEMON