1.1 --- a/lemon/matching.h Fri Sep 25 21:51:36 2009 +0200
1.2 +++ b/lemon/matching.h Sat Sep 26 10:17:31 2009 +0200
1.3 @@ -28,6 +28,7 @@
1.4 #include <lemon/unionfind.h>
1.5 #include <lemon/bin_heap.h>
1.6 #include <lemon/maps.h>
1.7 +#include <lemon/fractional_matching.h>
1.8
1.9 ///\ingroup matching
1.10 ///\file
1.11 @@ -797,6 +798,10 @@
1.12 BinHeap<Value, IntIntMap> *_delta4;
1.13
1.14 Value _delta_sum;
1.15 + int _unmatched;
1.16 +
1.17 + typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
1.18 + FractionalMatching *_fractional;
1.19
1.20 void createStructures() {
1.21 _node_num = countNodes(_graph);
1.22 @@ -1559,10 +1564,16 @@
1.23 _delta3_index(0), _delta3(0),
1.24 _delta4_index(0), _delta4(0),
1.25
1.26 - _delta_sum() {}
1.27 + _delta_sum(), _unmatched(0),
1.28 +
1.29 + _fractional(0)
1.30 + {}
1.31
1.32 ~MaxWeightedMatching() {
1.33 destroyStructures();
1.34 + if (_fractional) {
1.35 + delete _fractional;
1.36 + }
1.37 }
1.38
1.39 /// \name Execution Control
1.40 @@ -1591,6 +1602,8 @@
1.41 (*_delta4_index)[i] = _delta4->PRE_HEAP;
1.42 }
1.43
1.44 + _unmatched = _node_num;
1.45 +
1.46 int index = 0;
1.47 for (NodeIt n(_graph); n != INVALID; ++n) {
1.48 Value max = 0;
1.49 @@ -1625,18 +1638,155 @@
1.50 }
1.51 }
1.52
1.53 + /// \brief Initialize the algorithm with fractional matching
1.54 + ///
1.55 + /// This function initializes the algorithm with a fractional
1.56 + /// matching. This initialization is also called jumpstart heuristic.
1.57 + void fractionalInit() {
1.58 + createStructures();
1.59 +
1.60 + if (_fractional == 0) {
1.61 + _fractional = new FractionalMatching(_graph, _weight, false);
1.62 + }
1.63 + _fractional->run();
1.64 +
1.65 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.66 + (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1.67 + }
1.68 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.69 + (*_delta1_index)[n] = _delta1->PRE_HEAP;
1.70 + }
1.71 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.72 + (*_delta3_index)[e] = _delta3->PRE_HEAP;
1.73 + }
1.74 + for (int i = 0; i < _blossom_num; ++i) {
1.75 + (*_delta2_index)[i] = _delta2->PRE_HEAP;
1.76 + (*_delta4_index)[i] = _delta4->PRE_HEAP;
1.77 + }
1.78 +
1.79 + _unmatched = 0;
1.80 +
1.81 + int index = 0;
1.82 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.83 + Value pot = _fractional->nodeValue(n);
1.84 + (*_node_index)[n] = index;
1.85 + (*_node_data)[index].pot = pot;
1.86 + int blossom =
1.87 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
1.88 +
1.89 + (*_blossom_data)[blossom].status = MATCHED;
1.90 + (*_blossom_data)[blossom].pred = INVALID;
1.91 + (*_blossom_data)[blossom].next = _fractional->matching(n);
1.92 + if (_fractional->matching(n) == INVALID) {
1.93 + (*_blossom_data)[blossom].base = n;
1.94 + }
1.95 + (*_blossom_data)[blossom].pot = 0;
1.96 + (*_blossom_data)[blossom].offset = 0;
1.97 + ++index;
1.98 + }
1.99 +
1.100 + typename Graph::template NodeMap<bool> processed(_graph, false);
1.101 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.102 + if (processed[n]) continue;
1.103 + processed[n] = true;
1.104 + if (_fractional->matching(n) == INVALID) continue;
1.105 + int num = 1;
1.106 + Node v = _graph.target(_fractional->matching(n));
1.107 + while (n != v) {
1.108 + processed[v] = true;
1.109 + v = _graph.target(_fractional->matching(v));
1.110 + ++num;
1.111 + }
1.112 +
1.113 + if (num % 2 == 1) {
1.114 + std::vector<int> subblossoms(num);
1.115 +
1.116 + subblossoms[--num] = _blossom_set->find(n);
1.117 + _delta1->push(n, _fractional->nodeValue(n));
1.118 + v = _graph.target(_fractional->matching(n));
1.119 + while (n != v) {
1.120 + subblossoms[--num] = _blossom_set->find(v);
1.121 + _delta1->push(v, _fractional->nodeValue(v));
1.122 + v = _graph.target(_fractional->matching(v));
1.123 + }
1.124 +
1.125 + int surface =
1.126 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
1.127 + (*_blossom_data)[surface].status = EVEN;
1.128 + (*_blossom_data)[surface].pred = INVALID;
1.129 + (*_blossom_data)[surface].next = INVALID;
1.130 + (*_blossom_data)[surface].pot = 0;
1.131 + (*_blossom_data)[surface].offset = 0;
1.132 +
1.133 + _tree_set->insert(surface);
1.134 + ++_unmatched;
1.135 + }
1.136 + }
1.137 +
1.138 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.139 + int si = (*_node_index)[_graph.u(e)];
1.140 + int sb = _blossom_set->find(_graph.u(e));
1.141 + int ti = (*_node_index)[_graph.v(e)];
1.142 + int tb = _blossom_set->find(_graph.v(e));
1.143 + if ((*_blossom_data)[sb].status == EVEN &&
1.144 + (*_blossom_data)[tb].status == EVEN && sb != tb) {
1.145 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1.146 + dualScale * _weight[e]) / 2);
1.147 + }
1.148 + }
1.149 +
1.150 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.151 + int nb = _blossom_set->find(n);
1.152 + if ((*_blossom_data)[nb].status != MATCHED) continue;
1.153 + int ni = (*_node_index)[n];
1.154 +
1.155 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.156 + Node v = _graph.target(e);
1.157 + int vb = _blossom_set->find(v);
1.158 + int vi = (*_node_index)[v];
1.159 +
1.160 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.161 + dualScale * _weight[e];
1.162 +
1.163 + if ((*_blossom_data)[vb].status == EVEN) {
1.164 +
1.165 + int vt = _tree_set->find(vb);
1.166 +
1.167 + typename std::map<int, Arc>::iterator it =
1.168 + (*_node_data)[ni].heap_index.find(vt);
1.169 +
1.170 + if (it != (*_node_data)[ni].heap_index.end()) {
1.171 + if ((*_node_data)[ni].heap[it->second] > rw) {
1.172 + (*_node_data)[ni].heap.replace(it->second, e);
1.173 + (*_node_data)[ni].heap.decrease(e, rw);
1.174 + it->second = e;
1.175 + }
1.176 + } else {
1.177 + (*_node_data)[ni].heap.push(e, rw);
1.178 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
1.179 + }
1.180 + }
1.181 + }
1.182 +
1.183 + if (!(*_node_data)[ni].heap.empty()) {
1.184 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.185 + _delta2->push(nb, _blossom_set->classPrio(nb));
1.186 + }
1.187 + }
1.188 + }
1.189 +
1.190 /// \brief Start the algorithm
1.191 ///
1.192 /// This function starts the algorithm.
1.193 ///
1.194 - /// \pre \ref init() must be called before using this function.
1.195 + /// \pre \ref init() or \ref fractionalInit() must be called
1.196 + /// before using this function.
1.197 void start() {
1.198 enum OpType {
1.199 D1, D2, D3, D4
1.200 };
1.201
1.202 - int unmatched = _node_num;
1.203 - while (unmatched > 0) {
1.204 + while (_unmatched > 0) {
1.205 Value d1 = !_delta1->empty() ?
1.206 _delta1->prio() : std::numeric_limits<Value>::max();
1.207
1.208 @@ -1659,7 +1809,7 @@
1.209 {
1.210 Node n = _delta1->top();
1.211 unmatchNode(n);
1.212 - --unmatched;
1.213 + --_unmatched;
1.214 }
1.215 break;
1.216 case D2:
1.217 @@ -1669,7 +1819,7 @@
1.218 Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
1.219 if ((*_blossom_data)[blossom].next == INVALID) {
1.220 augmentOnArc(a);
1.221 - --unmatched;
1.222 + --_unmatched;
1.223 } else {
1.224 extendOnArc(a);
1.225 }
1.226 @@ -1692,7 +1842,7 @@
1.227 shrinkOnEdge(e, left_tree);
1.228 } else {
1.229 augmentOnEdge(e);
1.230 - unmatched -= 2;
1.231 + _unmatched -= 2;
1.232 }
1.233 }
1.234 } break;
1.235 @@ -1710,11 +1860,11 @@
1.236 ///
1.237 /// \note mwm.run() is just a shortcut of the following code.
1.238 /// \code
1.239 - /// mwm.init();
1.240 + /// mwm.fractionalInit();
1.241 /// mwm.start();
1.242 /// \endcode
1.243 void run() {
1.244 - init();
1.245 + fractionalInit();
1.246 start();
1.247 }
1.248
1.249 @@ -2074,6 +2224,11 @@
1.250 BinHeap<Value, IntIntMap> *_delta4;
1.251
1.252 Value _delta_sum;
1.253 + int _unmatched;
1.254 +
1.255 + typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
1.256 + FractionalMatching;
1.257 + FractionalMatching *_fractional;
1.258
1.259 void createStructures() {
1.260 _node_num = countNodes(_graph);
1.261 @@ -2789,10 +2944,16 @@
1.262 _delta3_index(0), _delta3(0),
1.263 _delta4_index(0), _delta4(0),
1.264
1.265 - _delta_sum() {}
1.266 + _delta_sum(), _unmatched(0),
1.267 +
1.268 + _fractional(0)
1.269 + {}
1.270
1.271 ~MaxWeightedPerfectMatching() {
1.272 destroyStructures();
1.273 + if (_fractional) {
1.274 + delete _fractional;
1.275 + }
1.276 }
1.277
1.278 /// \name Execution Control
1.279 @@ -2818,6 +2979,8 @@
1.280 (*_delta4_index)[i] = _delta4->PRE_HEAP;
1.281 }
1.282
1.283 + _unmatched = _node_num;
1.284 +
1.285 int index = 0;
1.286 for (NodeIt n(_graph); n != INVALID; ++n) {
1.287 Value max = - std::numeric_limits<Value>::max();
1.288 @@ -2851,18 +3014,152 @@
1.289 }
1.290 }
1.291
1.292 + /// \brief Initialize the algorithm with fractional matching
1.293 + ///
1.294 + /// This function initializes the algorithm with a fractional
1.295 + /// matching. This initialization is also called jumpstart heuristic.
1.296 + void fractionalInit() {
1.297 + createStructures();
1.298 +
1.299 + if (_fractional == 0) {
1.300 + _fractional = new FractionalMatching(_graph, _weight, false);
1.301 + }
1.302 + if (!_fractional->run()) {
1.303 + _unmatched = -1;
1.304 + return;
1.305 + }
1.306 +
1.307 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.308 + (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1.309 + }
1.310 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.311 + (*_delta3_index)[e] = _delta3->PRE_HEAP;
1.312 + }
1.313 + for (int i = 0; i < _blossom_num; ++i) {
1.314 + (*_delta2_index)[i] = _delta2->PRE_HEAP;
1.315 + (*_delta4_index)[i] = _delta4->PRE_HEAP;
1.316 + }
1.317 +
1.318 + _unmatched = 0;
1.319 +
1.320 + int index = 0;
1.321 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.322 + Value pot = _fractional->nodeValue(n);
1.323 + (*_node_index)[n] = index;
1.324 + (*_node_data)[index].pot = pot;
1.325 + int blossom =
1.326 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
1.327 +
1.328 + (*_blossom_data)[blossom].status = MATCHED;
1.329 + (*_blossom_data)[blossom].pred = INVALID;
1.330 + (*_blossom_data)[blossom].next = _fractional->matching(n);
1.331 + (*_blossom_data)[blossom].pot = 0;
1.332 + (*_blossom_data)[blossom].offset = 0;
1.333 + ++index;
1.334 + }
1.335 +
1.336 + typename Graph::template NodeMap<bool> processed(_graph, false);
1.337 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.338 + if (processed[n]) continue;
1.339 + processed[n] = true;
1.340 + if (_fractional->matching(n) == INVALID) continue;
1.341 + int num = 1;
1.342 + Node v = _graph.target(_fractional->matching(n));
1.343 + while (n != v) {
1.344 + processed[v] = true;
1.345 + v = _graph.target(_fractional->matching(v));
1.346 + ++num;
1.347 + }
1.348 +
1.349 + if (num % 2 == 1) {
1.350 + std::vector<int> subblossoms(num);
1.351 +
1.352 + subblossoms[--num] = _blossom_set->find(n);
1.353 + v = _graph.target(_fractional->matching(n));
1.354 + while (n != v) {
1.355 + subblossoms[--num] = _blossom_set->find(v);
1.356 + v = _graph.target(_fractional->matching(v));
1.357 + }
1.358 +
1.359 + int surface =
1.360 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
1.361 + (*_blossom_data)[surface].status = EVEN;
1.362 + (*_blossom_data)[surface].pred = INVALID;
1.363 + (*_blossom_data)[surface].next = INVALID;
1.364 + (*_blossom_data)[surface].pot = 0;
1.365 + (*_blossom_data)[surface].offset = 0;
1.366 +
1.367 + _tree_set->insert(surface);
1.368 + ++_unmatched;
1.369 + }
1.370 + }
1.371 +
1.372 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.373 + int si = (*_node_index)[_graph.u(e)];
1.374 + int sb = _blossom_set->find(_graph.u(e));
1.375 + int ti = (*_node_index)[_graph.v(e)];
1.376 + int tb = _blossom_set->find(_graph.v(e));
1.377 + if ((*_blossom_data)[sb].status == EVEN &&
1.378 + (*_blossom_data)[tb].status == EVEN && sb != tb) {
1.379 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1.380 + dualScale * _weight[e]) / 2);
1.381 + }
1.382 + }
1.383 +
1.384 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.385 + int nb = _blossom_set->find(n);
1.386 + if ((*_blossom_data)[nb].status != MATCHED) continue;
1.387 + int ni = (*_node_index)[n];
1.388 +
1.389 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.390 + Node v = _graph.target(e);
1.391 + int vb = _blossom_set->find(v);
1.392 + int vi = (*_node_index)[v];
1.393 +
1.394 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.395 + dualScale * _weight[e];
1.396 +
1.397 + if ((*_blossom_data)[vb].status == EVEN) {
1.398 +
1.399 + int vt = _tree_set->find(vb);
1.400 +
1.401 + typename std::map<int, Arc>::iterator it =
1.402 + (*_node_data)[ni].heap_index.find(vt);
1.403 +
1.404 + if (it != (*_node_data)[ni].heap_index.end()) {
1.405 + if ((*_node_data)[ni].heap[it->second] > rw) {
1.406 + (*_node_data)[ni].heap.replace(it->second, e);
1.407 + (*_node_data)[ni].heap.decrease(e, rw);
1.408 + it->second = e;
1.409 + }
1.410 + } else {
1.411 + (*_node_data)[ni].heap.push(e, rw);
1.412 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
1.413 + }
1.414 + }
1.415 + }
1.416 +
1.417 + if (!(*_node_data)[ni].heap.empty()) {
1.418 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.419 + _delta2->push(nb, _blossom_set->classPrio(nb));
1.420 + }
1.421 + }
1.422 + }
1.423 +
1.424 /// \brief Start the algorithm
1.425 ///
1.426 /// This function starts the algorithm.
1.427 ///
1.428 - /// \pre \ref init() must be called before using this function.
1.429 + /// \pre \ref init() or \ref fractionalInit() must be called before
1.430 + /// using this function.
1.431 bool start() {
1.432 enum OpType {
1.433 D2, D3, D4
1.434 };
1.435
1.436 - int unmatched = _node_num;
1.437 - while (unmatched > 0) {
1.438 + if (_unmatched == -1) return false;
1.439 +
1.440 + while (_unmatched > 0) {
1.441 Value d2 = !_delta2->empty() ?
1.442 _delta2->prio() : std::numeric_limits<Value>::max();
1.443
1.444 @@ -2906,7 +3203,7 @@
1.445 shrinkOnEdge(e, left_tree);
1.446 } else {
1.447 augmentOnEdge(e);
1.448 - unmatched -= 2;
1.449 + _unmatched -= 2;
1.450 }
1.451 }
1.452 } break;
1.453 @@ -2925,11 +3222,11 @@
1.454 ///
1.455 /// \note mwpm.run() is just a shortcut of the following code.
1.456 /// \code
1.457 - /// mwpm.init();
1.458 + /// mwpm.fractionalInit();
1.459 /// mwpm.start();
1.460 /// \endcode
1.461 bool run() {
1.462 - init();
1.463 + fractionalInit();
1.464 return start();
1.465 }
1.466
2.1 --- a/test/matching_test.cc Fri Sep 25 21:51:36 2009 +0200
2.2 +++ b/test/matching_test.cc Sat Sep 26 10:17:31 2009 +0200
2.3 @@ -401,22 +401,46 @@
2.4 graphReader(graph, lgfs).
2.5 edgeMap("weight", weight).run();
2.6
2.7 - MaxMatching<SmartGraph> mm(graph);
2.8 - mm.run();
2.9 - checkMatching(graph, mm);
2.10 + bool perfect;
2.11 + {
2.12 + MaxMatching<SmartGraph> mm(graph);
2.13 + mm.run();
2.14 + checkMatching(graph, mm);
2.15 + perfect = 2 * mm.matchingSize() == countNodes(graph);
2.16 + }
2.17
2.18 - MaxWeightedMatching<SmartGraph> mwm(graph, weight);
2.19 - mwm.run();
2.20 - checkWeightedMatching(graph, weight, mwm);
2.21 + {
2.22 + MaxWeightedMatching<SmartGraph> mwm(graph, weight);
2.23 + mwm.run();
2.24 + checkWeightedMatching(graph, weight, mwm);
2.25 + }
2.26
2.27 - MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
2.28 - bool perfect = mwpm.run();
2.29 + {
2.30 + MaxWeightedMatching<SmartGraph> mwm(graph, weight);
2.31 + mwm.init();
2.32 + mwm.start();
2.33 + checkWeightedMatching(graph, weight, mwm);
2.34 + }
2.35
2.36 - check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
2.37 - "Perfect matching found");
2.38 + {
2.39 + MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
2.40 + bool result = mwpm.run();
2.41 +
2.42 + check(result == perfect, "Perfect matching found");
2.43 + if (perfect) {
2.44 + checkWeightedPerfectMatching(graph, weight, mwpm);
2.45 + }
2.46 + }
2.47
2.48 - if (perfect) {
2.49 - checkWeightedPerfectMatching(graph, weight, mwpm);
2.50 + {
2.51 + MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
2.52 + mwpm.init();
2.53 + bool result = mwpm.start();
2.54 +
2.55 + check(result == perfect, "Perfect matching found");
2.56 + if (perfect) {
2.57 + checkWeightedPerfectMatching(graph, weight, mwpm);
2.58 + }
2.59 }
2.60 }
2.61