lemon: comparison lemon/matching.h

equal deleted inserted replaced

-:4fa2e7cacf26
+:b397b2f8308e
 #include <lemon/core.h>
 #include <lemon/unionfind.h>
 #include <lemon/bin_heap.h>
 #include <lemon/maps.h>
+#include <lemon/fractional_matching.h>
 ///\ingroup matching
 ///\file
 ///\brief Maximum matching algorithms in general graphs.
 IntIntMap *_delta4_index;
 BinHeap<Value, IntIntMap> *_delta4;
 Value _delta_sum;
+int _unmatched;
+typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
+FractionalMatching *_fractional;
 void createStructures() {
 _node_num = countNodes(_graph);
 _blossom_num = _node_num * 3 / 2;
 _delta1_index(0), _delta1(0),
 _delta2_index(0), _delta2(0),
 _delta3_index(0), _delta3(0),
 _delta4_index(0), _delta4(0),
-_delta_sum() {}
+_delta_sum(), _unmatched(0),
+_fractional(0)
+{}
 ~MaxWeightedMatching() {
 destroyStructures();
+if (_fractional) {
+delete _fractional;
+}
 }
 /// \name Execution Control
 /// The simplest way to execute the algorithm is to use the
 /// \ref run() member function.
 }
 for (int i = 0; i < _blossom_num; ++i) {
 (*_delta2_index)[i] = _delta2->PRE_HEAP;
 (*_delta4_index)[i] = _delta4->PRE_HEAP;
 }
+_unmatched = _node_num;
 int index = 0;
 for (NodeIt n(_graph); n != INVALID; ++n) {
 Value max = 0;
 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 dualScale * _weight[e]) / 2);
 }
 }
 }
+/// \brief Initialize the algorithm with fractional matching
+///
+/// This function initializes the algorithm with a fractional
+/// matching. This initialization is also called jumpstart heuristic.
+void fractionalInit() {
+createStructures();
+if (_fractional == 0) {
+_fractional = new FractionalMatching(_graph, _weight, false);
+}
+_fractional->run();
+for (ArcIt e(_graph); e != INVALID; ++e) {
+(*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
+}
+for (NodeIt n(_graph); n != INVALID; ++n) {
+(*_delta1_index)[n] = _delta1->PRE_HEAP;
+}
+for (EdgeIt e(_graph); e != INVALID; ++e) {
+(*_delta3_index)[e] = _delta3->PRE_HEAP;
+}
+for (int i = 0; i < _blossom_num; ++i) {
+(*_delta2_index)[i] = _delta2->PRE_HEAP;
+(*_delta4_index)[i] = _delta4->PRE_HEAP;
+}
+_unmatched = 0;
+int index = 0;
+for (NodeIt n(_graph); n != INVALID; ++n) {
+Value pot = _fractional->nodeValue(n);
+(*_node_index)[n] = index;
+(*_node_data)[index].pot = pot;
+int blossom =
+_blossom_set->insert(n, std::numeric_limits<Value>::max());
+(*_blossom_data)[blossom].status = MATCHED;
+(*_blossom_data)[blossom].pred = INVALID;
+(*_blossom_data)[blossom].next = _fractional->matching(n);
+if (_fractional->matching(n) == INVALID) {
+(*_blossom_data)[blossom].base = n;
+}
+(*_blossom_data)[blossom].pot = 0;
+(*_blossom_data)[blossom].offset = 0;
+++index;
+}
+typename Graph::template NodeMap<bool> processed(_graph, false);
+for (NodeIt n(_graph); n != INVALID; ++n) {
+if (processed[n]) continue;
+processed[n] = true;
+if (_fractional->matching(n) == INVALID) continue;
+int num = 1;
+Node v = _graph.target(_fractional->matching(n));
+while (n != v) {
+processed[v] = true;
+v = _graph.target(_fractional->matching(v));
+++num;
+}
+if (num % 2 == 1) {
+std::vector<int> subblossoms(num);
+subblossoms[--num] = _blossom_set->find(n);
+_delta1->push(n, _fractional->nodeValue(n));
+v = _graph.target(_fractional->matching(n));
+while (n != v) {
+subblossoms[--num] = _blossom_set->find(v);
+_delta1->push(v, _fractional->nodeValue(v));
+v = _graph.target(_fractional->matching(v));
+}
+int surface =
+_blossom_set->join(subblossoms.begin(), subblossoms.end());
+(*_blossom_data)[surface].status = EVEN;
+(*_blossom_data)[surface].pred = INVALID;
+(*_blossom_data)[surface].next = INVALID;
+(*_blossom_data)[surface].pot = 0;
+(*_blossom_data)[surface].offset = 0;
+_tree_set->insert(surface);
+++_unmatched;
+}
+}
+for (EdgeIt e(_graph); e != INVALID; ++e) {
+int si = (*_node_index)[_graph.u(e)];
+int sb = _blossom_set->find(_graph.u(e));
+int ti = (*_node_index)[_graph.v(e)];
+int tb = _blossom_set->find(_graph.v(e));
+if ((*_blossom_data)[sb].status == EVEN &&
+(*_blossom_data)[tb].status == EVEN && sb != tb) {
+_delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
+dualScale * _weight[e]) / 2);
+}
+}
+for (NodeIt n(_graph); n != INVALID; ++n) {
+int nb = _blossom_set->find(n);
+if ((*_blossom_data)[nb].status != MATCHED) continue;
+int ni = (*_node_index)[n];
+for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+Node v = _graph.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) {
+int vt = _tree_set->find(vb);
+typename std::map<int, Arc>::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, e);
+(*_node_data)[ni].heap.decrease(e, rw);
+it->second = e;
+}
+} else {
+(*_node_data)[ni].heap.push(e, rw);
+(*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
+}
+}
+}
+if (!(*_node_data)[ni].heap.empty()) {
+_blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
+_delta2->push(nb, _blossom_set->classPrio(nb));
+}
+}
+}
 /// \brief Start the algorithm
 ///
 /// This function starts the algorithm.
 ///
-/// \pre \ref init() must be called before using this function.
+/// \pre \ref init() or \ref fractionalInit() must be called
+/// before using this function.
 void start() {
 enum OpType {
 D1, D2, D3, D4
 };
-int unmatched = _node_num;
+while (_unmatched > 0) {
-while (unmatched > 0) {
 Value d1 = !_delta1->empty() ?
 _delta1->prio() : std::numeric_limits<Value>::max();
 Value d2 = !_delta2->empty() ?
 _delta2->prio() : std::numeric_limits<Value>::max();
 switch (ot) {
 case D1:
 {
 Node n = _delta1->top();
 unmatchNode(n);
---unmatched;
+--_unmatched;
 }
 break;
 case D2:
 {
 int blossom = _delta2->top();
 Node n = _blossom_set->classTop(blossom);
 Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
 if ((*_blossom_data)[blossom].next == INVALID) {
 augmentOnArc(a);
---unmatched;
+--_unmatched;
 } else {
 extendOnArc(a);
 }
 }
 break;
 if (left_tree == right_tree) {
 shrinkOnEdge(e, left_tree);
 } else {
 augmentOnEdge(e);
-unmatched -= 2;
+_unmatched -= 2;
 }
 }
 } break;
 case D4:
 splitBlossom(_delta4->top());
 ///
 /// This method runs the \c %MaxWeightedMatching algorithm.
 ///
 /// \note mwm.run() is just a shortcut of the following code.
 /// \code
-///   mwm.init();
+///   mwm.fractionalInit();
 ///   mwm.start();
 /// \endcode
 void run() {
-init();
+fractionalInit();
 start();
 }
 /// @}
 IntIntMap *_delta4_index;
 BinHeap<Value, IntIntMap> *_delta4;
 Value _delta_sum;
+int _unmatched;
+typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
+FractionalMatching;
+FractionalMatching *_fractional;
 void createStructures() {
 _node_num = countNodes(_graph);
 _blossom_num = _node_num * 3 / 2;
 _delta2_index(0), _delta2(0),
 _delta3_index(0), _delta3(0),
 _delta4_index(0), _delta4(0),
-_delta_sum() {}
+_delta_sum(), _unmatched(0),
+_fractional(0)
+{}
 ~MaxWeightedPerfectMatching() {
 destroyStructures();
+if (_fractional) {
+delete _fractional;
+}
 }
 /// \name Execution Control
 /// The simplest way to execute the algorithm is to use the
 /// \ref run() member function.
 }
 for (int i = 0; i < _blossom_num; ++i) {
 (*_delta2_index)[i] = _delta2->PRE_HEAP;
 (*_delta4_index)[i] = _delta4->PRE_HEAP;
 }
+_unmatched = _node_num;
 int index = 0;
 for (NodeIt n(_graph); n != INVALID; ++n) {
 Value max = - std::numeric_limits<Value>::max();
 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 dualScale * _weight[e]) / 2);
 }
 }
 }
+/// \brief Initialize the algorithm with fractional matching
+///
+/// This function initializes the algorithm with a fractional
+/// matching. This initialization is also called jumpstart heuristic.
+void fractionalInit() {
+createStructures();
+if (_fractional == 0) {
+_fractional = new FractionalMatching(_graph, _weight, false);
+}
+if (!_fractional->run()) {
+_unmatched = -1;
+return;
+}
+for (ArcIt e(_graph); e != INVALID; ++e) {
+(*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
+}
+for (EdgeIt e(_graph); e != INVALID; ++e) {
+(*_delta3_index)[e] = _delta3->PRE_HEAP;
+}
+for (int i = 0; i < _blossom_num; ++i) {
+(*_delta2_index)[i] = _delta2->PRE_HEAP;
+(*_delta4_index)[i] = _delta4->PRE_HEAP;
+}
+_unmatched = 0;
+int index = 0;
+for (NodeIt n(_graph); n != INVALID; ++n) {
+Value pot = _fractional->nodeValue(n);
+(*_node_index)[n] = index;
+(*_node_data)[index].pot = pot;
+int blossom =
+_blossom_set->insert(n, std::numeric_limits<Value>::max());
+(*_blossom_data)[blossom].status = MATCHED;
+(*_blossom_data)[blossom].pred = INVALID;
+(*_blossom_data)[blossom].next = _fractional->matching(n);
+(*_blossom_data)[blossom].pot = 0;
+(*_blossom_data)[blossom].offset = 0;
+++index;
+}
+typename Graph::template NodeMap<bool> processed(_graph, false);
+for (NodeIt n(_graph); n != INVALID; ++n) {
+if (processed[n]) continue;
+processed[n] = true;
+if (_fractional->matching(n) == INVALID) continue;
+int num = 1;
+Node v = _graph.target(_fractional->matching(n));
+while (n != v) {
+processed[v] = true;
+v = _graph.target(_fractional->matching(v));
+++num;
+}
+if (num % 2 == 1) {
+std::vector<int> subblossoms(num);
+subblossoms[--num] = _blossom_set->find(n);
+v = _graph.target(_fractional->matching(n));
+while (n != v) {
+subblossoms[--num] = _blossom_set->find(v);
+v = _graph.target(_fractional->matching(v));
+}
+int surface =
+_blossom_set->join(subblossoms.begin(), subblossoms.end());
+(*_blossom_data)[surface].status = EVEN;
+(*_blossom_data)[surface].pred = INVALID;
+(*_blossom_data)[surface].next = INVALID;
+(*_blossom_data)[surface].pot = 0;
+(*_blossom_data)[surface].offset = 0;
+_tree_set->insert(surface);
+++_unmatched;
+}
+}
+for (EdgeIt e(_graph); e != INVALID; ++e) {
+int si = (*_node_index)[_graph.u(e)];
+int sb = _blossom_set->find(_graph.u(e));
+int ti = (*_node_index)[_graph.v(e)];
+int tb = _blossom_set->find(_graph.v(e));
+if ((*_blossom_data)[sb].status == EVEN &&
+(*_blossom_data)[tb].status == EVEN && sb != tb) {
+_delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
+dualScale * _weight[e]) / 2);
+}
+}
+for (NodeIt n(_graph); n != INVALID; ++n) {
+int nb = _blossom_set->find(n);
+if ((*_blossom_data)[nb].status != MATCHED) continue;
+int ni = (*_node_index)[n];
+for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+Node v = _graph.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) {
+int vt = _tree_set->find(vb);
+typename std::map<int, Arc>::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, e);
+(*_node_data)[ni].heap.decrease(e, rw);
+it->second = e;
+}
+} else {
+(*_node_data)[ni].heap.push(e, rw);
+(*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
+}
+}
+}
+if (!(*_node_data)[ni].heap.empty()) {
+_blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
+_delta2->push(nb, _blossom_set->classPrio(nb));
+}
+}
+}
 /// \brief Start the algorithm
 ///
 /// This function starts the algorithm.
 ///
-/// \pre \ref init() must be called before using this function.
+/// \pre \ref init() or \ref fractionalInit() must be called before
+/// using this function.
 bool start() {
 enum OpType {
 D2, D3, D4
 };
-int unmatched = _node_num;
+if (_unmatched == -1) return false;
-while (unmatched > 0) {
+while (_unmatched > 0) {
 Value d2 = !_delta2->empty() ?
 _delta2->prio() : std::numeric_limits<Value>::max();
 Value d3 = !_delta3->empty() ?
 _delta3->prio() : std::numeric_limits<Value>::max();
 if (left_tree == right_tree) {
 shrinkOnEdge(e, left_tree);
 } else {
 augmentOnEdge(e);
-unmatched -= 2;
+_unmatched -= 2;
 }
 }
 } break;
 case D4:
 splitBlossom(_delta4->top());
 ///
 /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
 ///
 /// \note mwpm.run() is just a shortcut of the following code.
 /// \code
-///   mwpm.init();
+///   mwpm.fractionalInit();
 ///   mwpm.start();
 /// \endcode
 bool run() {
-init();
+fractionalInit();
 return start();
 }
 /// @}

changeset 953	d8ea85825e02
parent 947	0513ccfea967
child 954	07ec2b52e53d