gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Imporvements for the matching algorithms (#264)
0 3 0
default
3 files changed with 393 insertions and 188 deletions:
↑ Collapse diff ↑
Ignore white space 6 line context
... ...
@@ -246,387 +246,388 @@
246 246
/**
247 247
@defgroup paths Path Structures
248 248
@ingroup datas
249 249
\brief %Path structures implemented in LEMON.
250 250

	
251 251
This group contains the path structures implemented in LEMON.
252 252

	
253 253
LEMON provides flexible data structures to work with paths.
254 254
All of them have similar interfaces and they can be copied easily with
255 255
assignment operators and copy constructors. This makes it easy and
256 256
efficient to have e.g. the Dijkstra algorithm to store its result in
257 257
any kind of path structure.
258 258

	
259 259
\sa lemon::concepts::Path
260 260
*/
261 261

	
262 262
/**
263 263
@defgroup auxdat Auxiliary Data Structures
264 264
@ingroup datas
265 265
\brief Auxiliary data structures implemented in LEMON.
266 266

	
267 267
This group contains some data structures implemented in LEMON in
268 268
order to make it easier to implement combinatorial algorithms.
269 269
*/
270 270

	
271 271
/**
272 272
@defgroup algs Algorithms
273 273
\brief This group contains the several algorithms
274 274
implemented in LEMON.
275 275

	
276 276
This group contains the several algorithms
277 277
implemented in LEMON.
278 278
*/
279 279

	
280 280
/**
281 281
@defgroup search Graph Search
282 282
@ingroup algs
283 283
\brief Common graph search algorithms.
284 284

	
285 285
This group contains the common graph search algorithms, namely
286 286
\e breadth-first \e search (BFS) and \e depth-first \e search (DFS).
287 287
*/
288 288

	
289 289
/**
290 290
@defgroup shortest_path Shortest Path Algorithms
291 291
@ingroup algs
292 292
\brief Algorithms for finding shortest paths.
293 293

	
294 294
This group contains the algorithms for finding shortest paths in digraphs.
295 295

	
296 296
 - \ref Dijkstra algorithm for finding shortest paths from a source node
297 297
   when all arc lengths are non-negative.
298 298
 - \ref BellmanFord "Bellman-Ford" algorithm for finding shortest paths
299 299
   from a source node when arc lenghts can be either positive or negative,
300 300
   but the digraph should not contain directed cycles with negative total
301 301
   length.
302 302
 - \ref FloydWarshall "Floyd-Warshall" and \ref Johnson "Johnson" algorithms
303 303
   for solving the \e all-pairs \e shortest \e paths \e problem when arc
304 304
   lenghts can be either positive or negative, but the digraph should
305 305
   not contain directed cycles with negative total length.
306 306
 - \ref Suurballe A successive shortest path algorithm for finding
307 307
   arc-disjoint paths between two nodes having minimum total length.
308 308
*/
309 309

	
310 310
/**
311 311
@defgroup max_flow Maximum Flow Algorithms
312 312
@ingroup algs
313 313
\brief Algorithms for finding maximum flows.
314 314

	
315 315
This group contains the algorithms for finding maximum flows and
316 316
feasible circulations.
317 317

	
318 318
The \e maximum \e flow \e problem is to find a flow of maximum value between
319 319
a single source and a single target. Formally, there is a \f$G=(V,A)\f$
320 320
digraph, a \f$cap:A\rightarrow\mathbf{R}^+_0\f$ capacity function and
321 321
\f$s, t \in V\f$ source and target nodes.
322 322
A maximum flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of the
323 323
following optimization problem.
324 324

	
325 325
\f[ \max\sum_{a\in\delta_{out}(s)}f(a) - \sum_{a\in\delta_{in}(s)}f(a) \f]
326 326
\f[ \sum_{a\in\delta_{out}(v)} f(a) = \sum_{a\in\delta_{in}(v)} f(a)
327 327
    \qquad \forall v\in V\setminus\{s,t\} \f]
328 328
\f[ 0 \leq f(a) \leq cap(a) \qquad \forall a\in A \f]
329 329

	
330 330
LEMON contains several algorithms for solving maximum flow problems:
331 331
- \ref EdmondsKarp Edmonds-Karp algorithm.
332 332
- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm.
333 333
- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees.
334 334
- \ref GoldbergTarjan Preflow push-relabel algorithm with dynamic trees.
335 335

	
336 336
In most cases the \ref Preflow "Preflow" algorithm provides the
337 337
fastest method for computing a maximum flow. All implementations
338 338
provides functions to also query the minimum cut, which is the dual
339 339
problem of the maximum flow.
340 340
*/
341 341

	
342 342
/**
343 343
@defgroup min_cost_flow Minimum Cost Flow Algorithms
344 344
@ingroup algs
345 345

	
346 346
\brief Algorithms for finding minimum cost flows and circulations.
347 347

	
348 348
This group contains the algorithms for finding minimum cost flows and
349 349
circulations.
350 350

	
351 351
The \e minimum \e cost \e flow \e problem is to find a feasible flow of
352 352
minimum total cost from a set of supply nodes to a set of demand nodes
353 353
in a network with capacity constraints and arc costs.
354 354
Formally, let \f$G=(V,A)\f$ be a digraph,
355 355
\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
356 356
upper bounds for the flow values on the arcs,
357 357
\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
358 358
on the arcs, and
359 359
\f$supply: V\rightarrow\mathbf{Z}\f$ denotes the supply/demand values
360 360
of the nodes.
361 361
A minimum cost flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of
362 362
the following optimization problem.
363 363

	
364 364
\f[ \min\sum_{a\in A} f(a) cost(a) \f]
365 365
\f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a) =
366 366
    supply(v) \qquad \forall v\in V \f]
367 367
\f[ lower(a) \leq f(a) \leq upper(a) \qquad \forall a\in A \f]
368 368

	
369 369
LEMON contains several algorithms for solving minimum cost flow problems:
370 370
 - \ref CycleCanceling Cycle-canceling algorithms.
371 371
 - \ref CapacityScaling Successive shortest path algorithm with optional
372 372
   capacity scaling.
373 373
 - \ref CostScaling Push-relabel and augment-relabel algorithms based on
374 374
   cost scaling.
375 375
 - \ref NetworkSimplex Primal network simplex algorithm with various
376 376
   pivot strategies.
377 377
*/
378 378

	
379 379
/**
380 380
@defgroup min_cut Minimum Cut Algorithms
381 381
@ingroup algs
382 382

	
383 383
\brief Algorithms for finding minimum cut in graphs.
384 384

	
385 385
This group contains the algorithms for finding minimum cut in graphs.
386 386

	
387 387
The \e minimum \e cut \e problem is to find a non-empty and non-complete
388 388
\f$X\f$ subset of the nodes with minimum overall capacity on
389 389
outgoing arcs. Formally, there is a \f$G=(V,A)\f$ digraph, a
390 390
\f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function. The minimum
391 391
cut is the \f$X\f$ solution of the next optimization problem:
392 392

	
393 393
\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}}
394 394
    \sum_{uv\in A, u\in X, v\not\in X}cap(uv) \f]
395 395

	
396 396
LEMON contains several algorithms related to minimum cut problems:
397 397

	
398 398
- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut
399 399
  in directed graphs.
400 400
- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for
401 401
  calculating minimum cut in undirected graphs.
402 402
- \ref GomoryHu "Gomory-Hu tree computation" for calculating
403 403
  all-pairs minimum cut in undirected graphs.
404 404

	
405 405
If you want to find minimum cut just between two distinict nodes,
406 406
see the \ref max_flow "maximum flow problem".
407 407
*/
408 408

	
409 409
/**
410 410
@defgroup graph_properties Connectivity and Other Graph Properties
411 411
@ingroup algs
412 412
\brief Algorithms for discovering the graph properties
413 413

	
414 414
This group contains the algorithms for discovering the graph properties
415 415
like connectivity, bipartiteness, euler property, simplicity etc.
416 416

	
417 417
\image html edge_biconnected_components.png
418 418
\image latex edge_biconnected_components.eps "bi-edge-connected components" width=\textwidth
419 419
*/
420 420

	
421 421
/**
422 422
@defgroup planar Planarity Embedding and Drawing
423 423
@ingroup algs
424 424
\brief Algorithms for planarity checking, embedding and drawing
425 425

	
426 426
This group contains the algorithms for planarity checking,
427 427
embedding and drawing.
428 428

	
429 429
\image html planar.png
430 430
\image latex planar.eps "Plane graph" width=\textwidth
431 431
*/
432 432

	
433 433
/**
434 434
@defgroup matching Matching Algorithms
435 435
@ingroup algs
436 436
\brief Algorithms for finding matchings in graphs and bipartite graphs.
437 437

	
438
This group contains algorithm objects and functions to calculate
438
This group contains the algorithms for calculating
439 439
matchings in graphs and bipartite graphs. The general matching problem is
440
finding a subset of the arcs which does not shares common endpoints.
440
finding a subset of the edges for which each node has at most one incident
441
edge.
441 442

	
442 443
There are several different algorithms for calculate matchings in
443 444
graphs.  The matching problems in bipartite graphs are generally
444 445
easier than in general graphs. The goal of the matching optimization
445 446
can be finding maximum cardinality, maximum weight or minimum cost
446 447
matching. The search can be constrained to find perfect or
447 448
maximum cardinality matching.
448 449

	
449 450
The matching algorithms implemented in LEMON:
450 451
- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm
451 452
  for calculating maximum cardinality matching in bipartite graphs.
452 453
- \ref PrBipartiteMatching Push-relabel algorithm
453 454
  for calculating maximum cardinality matching in bipartite graphs.
454 455
- \ref MaxWeightedBipartiteMatching
455 456
  Successive shortest path algorithm for calculating maximum weighted
456 457
  matching and maximum weighted bipartite matching in bipartite graphs.
457 458
- \ref MinCostMaxBipartiteMatching
458 459
  Successive shortest path algorithm for calculating minimum cost maximum
459 460
  matching in bipartite graphs.
460 461
- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating
461 462
  maximum cardinality matching in general graphs.
462 463
- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating
463 464
  maximum weighted matching in general graphs.
464 465
- \ref MaxWeightedPerfectMatching
465 466
  Edmond's blossom shrinking algorithm for calculating maximum weighted
466 467
  perfect matching in general graphs.
467 468

	
468 469
\image html bipartite_matching.png
469 470
\image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth
470 471
*/
471 472

	
472 473
/**
473 474
@defgroup spantree Minimum Spanning Tree Algorithms
474 475
@ingroup algs
475 476
\brief Algorithms for finding a minimum cost spanning tree in a graph.
476 477

	
477 478
This group contains the algorithms for finding a minimum cost spanning
478 479
tree in a graph.
479 480
*/
480 481

	
481 482
/**
482 483
@defgroup auxalg Auxiliary Algorithms
483 484
@ingroup algs
484 485
\brief Auxiliary algorithms implemented in LEMON.
485 486

	
486 487
This group contains some algorithms implemented in LEMON
487 488
in order to make it easier to implement complex algorithms.
488 489
*/
489 490

	
490 491
/**
491 492
@defgroup approx Approximation Algorithms
492 493
@ingroup algs
493 494
\brief Approximation algorithms.
494 495

	
495 496
This group contains the approximation and heuristic algorithms
496 497
implemented in LEMON.
497 498
*/
498 499

	
499 500
/**
500 501
@defgroup gen_opt_group General Optimization Tools
501 502
\brief This group contains some general optimization frameworks
502 503
implemented in LEMON.
503 504

	
504 505
This group contains some general optimization frameworks
505 506
implemented in LEMON.
506 507
*/
507 508

	
508 509
/**
509 510
@defgroup lp_group Lp and Mip Solvers
510 511
@ingroup gen_opt_group
511 512
\brief Lp and Mip solver interfaces for LEMON.
512 513

	
513 514
This group contains Lp and Mip solver interfaces for LEMON. The
514 515
various LP solvers could be used in the same manner with this
515 516
interface.
516 517
*/
517 518

	
518 519
/**
519 520
@defgroup lp_utils Tools for Lp and Mip Solvers
520 521
@ingroup lp_group
521 522
\brief Helper tools to the Lp and Mip solvers.
522 523

	
523 524
This group adds some helper tools to general optimization framework
524 525
implemented in LEMON.
525 526
*/
526 527

	
527 528
/**
528 529
@defgroup metah Metaheuristics
529 530
@ingroup gen_opt_group
530 531
\brief Metaheuristics for LEMON library.
531 532

	
532 533
This group contains some metaheuristic optimization tools.
533 534
*/
534 535

	
535 536
/**
536 537
@defgroup utils Tools and Utilities
537 538
\brief Tools and utilities for programming in LEMON
538 539

	
539 540
Tools and utilities for programming in LEMON.
540 541
*/
541 542

	
542 543
/**
543 544
@defgroup gutils Basic Graph Utilities
544 545
@ingroup utils
545 546
\brief Simple basic graph utilities.
546 547

	
547 548
This group contains some simple basic graph utilities.
548 549
*/
549 550

	
550 551
/**
551 552
@defgroup misc Miscellaneous Tools
552 553
@ingroup utils
553 554
\brief Tools for development, debugging and testing.
554 555

	
555 556
This group contains several useful tools for development,
556 557
debugging and testing.
557 558
*/
558 559

	
559 560
/**
560 561
@defgroup timecount Time Measuring and Counting
561 562
@ingroup misc
562 563
\brief Simple tools for measuring the performance of algorithms.
563 564

	
564 565
This group contains simple tools for measuring the performance
565 566
of algorithms.
566 567
*/
567 568

	
568 569
/**
569 570
@defgroup exceptions Exceptions
570 571
@ingroup utils
571 572
\brief Exceptions defined in LEMON.
572 573

	
573 574
This group contains the exceptions defined in LEMON.
574 575
*/
575 576

	
576 577
/**
577 578
@defgroup io_group Input-Output
578 579
\brief Graph Input-Output methods
579 580

	
580 581
This group contains the tools for importing and exporting graphs
581 582
and graph related data. Now it supports the \ref lgf-format
582 583
"LEMON Graph Format", the \c DIMACS format and the encapsulated
583 584
postscript (EPS) format.
584 585
*/
585 586

	
586 587
/**
587 588
@defgroup lemon_io LEMON Graph Format
588 589
@ingroup io_group
589 590
\brief Reading and writing LEMON Graph Format.
590 591

	
591 592
This group contains methods for reading and writing
592 593
\ref lgf-format "LEMON Graph Format".
593 594
*/
594 595

	
595 596
/**
596 597
@defgroup eps_io Postscript Exporting
597 598
@ingroup io_group
598 599
\brief General \c EPS drawer and graph exporter
599 600

	
600 601
This group contains general \c EPS drawing methods and special
601 602
graph exporting tools.
602 603
*/
603 604

	
604 605
/**
605 606
@defgroup dimacs_group DIMACS format
606 607
@ingroup io_group
607 608
\brief Read and write files in DIMACS format
608 609

	
609 610
Tools to read a digraph from or write it to a file in DIMACS format data.
610 611
*/
611 612

	
612 613
/**
613 614
@defgroup nauty_group NAUTY Format
614 615
@ingroup io_group
615 616
\brief Read \e Nauty format
616 617

	
617 618
Tool to read graphs from \e Nauty format data.
618 619
*/
619 620

	
620 621
/**
621 622
@defgroup concept Concepts
622 623
\brief Skeleton classes and concept checking classes
623 624

	
624 625
This group contains the data/algorithm skeletons and concept checking
625 626
classes implemented in LEMON.
626 627

	
627 628
The purpose of the classes in this group is fourfold.
628 629

	
629 630
- These classes contain the documentations of the %concepts. In order
630 631
  to avoid document multiplications, an implementation of a concept
631 632
  simply refers to the corresponding concept class.
632 633

	
Ignore white space 384 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_MAX_MATCHING_H
20 20
#define LEMON_MAX_MATCHING_H
21 21

	
22 22
#include <vector>
23 23
#include <queue>
24 24
#include <set>
25 25
#include <limits>
26 26

	
27 27
#include <lemon/core.h>
28 28
#include <lemon/unionfind.h>
29 29
#include <lemon/bin_heap.h>
30 30
#include <lemon/maps.h>
31 31

	
32 32
///\ingroup matching
33 33
///\file
34 34
///\brief Maximum matching algorithms in general graphs.
35 35

	
36 36
namespace lemon {
37 37

	
38 38
  /// \ingroup matching
39 39
  ///
40
  /// \brief Edmonds' alternating forest maximum matching algorithm.
40
  /// \brief Maximum cardinality matching in general graphs
41 41
  ///
42
  /// This class implements Edmonds' alternating forest matching
43
  /// algorithm. The algorithm can be started from an arbitrary initial
44
  /// matching (the default is the empty one)
42
  /// This class implements Edmonds' alternating forest matching algorithm
43
  /// for finding a maximum cardinality matching in a general graph. 
44
  /// It can be started from an arbitrary initial matching 
45
  /// (the default is the empty one).
45 46
  ///
46 47
  /// The dual solution of the problem is a map of the nodes to
47
  /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
48
  /// MATCHED/C showing the Gallai-Edmonds decomposition of the
49
  /// graph. The nodes in \c EVEN/D induce a graph with
50
  /// factor-critical components, the nodes in \c ODD/A form the
51
  /// barrier, and the nodes in \c MATCHED/C induce a graph having a
52
  /// perfect matching. The number of the factor-critical components
48
  /// \ref MaxMatching::Status "Status", having values \c EVEN (or \c D),
49
  /// \c ODD (or \c A) and \c MATCHED (or \c C) defining the Gallai-Edmonds
50
  /// decomposition of the graph. The nodes in \c EVEN/D induce a subgraph
51
  /// with factor-critical components, the nodes in \c ODD/A form the
52
  /// canonical barrier, and the nodes in \c MATCHED/C induce a graph having
53
  /// a perfect matching. The number of the factor-critical components
53 54
  /// minus the number of barrier nodes is a lower bound on the
54 55
  /// unmatched nodes, and the matching is optimal if and only if this bound is
55
  /// tight. This decomposition can be attained by calling \c
56
  /// tight. This decomposition can be obtained by calling \c
56 57
  /// decomposition() after running the algorithm.
57 58
  ///
58
  /// \param GR The graph type the algorithm runs on.
59
  /// \tparam GR The graph type the algorithm runs on.
59 60
  template <typename GR>
60 61
  class MaxMatching {
61 62
  public:
62 63

	
64
    /// The graph type of the algorithm
63 65
    typedef GR Graph;
64 66
    typedef typename Graph::template NodeMap<typename Graph::Arc>
65 67
    MatchingMap;
66 68

	
67
    ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
69
    ///\brief Status constants for Gallai-Edmonds decomposition.
68 70
    ///
69
    ///Indicates the Gallai-Edmonds decomposition of the graph. The
70
    ///nodes with Status \c EVEN/D induce a graph with factor-critical
71
    ///components, the nodes in \c ODD/A form the canonical barrier,
72
    ///and the nodes in \c MATCHED/C induce a graph having a perfect
73
    ///matching.
71
    ///These constants are used for indicating the Gallai-Edmonds 
72
    ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
73
    ///induce a subgraph with factor-critical components, the nodes with
74
    ///status \c ODD (or \c A) form the canonical barrier, and the nodes
75
    ///with status \c MATCHED (or \c C) induce a subgraph having a 
76
    ///perfect matching.
74 77
    enum Status {
75
      EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
78
      EVEN = 1,       ///< = 1. (\c D is an alias for \c EVEN.)
79
      D = 1,
80
      MATCHED = 0,    ///< = 0. (\c C is an alias for \c MATCHED.)
81
      C = 0,
82
      ODD = -1,       ///< = -1. (\c A is an alias for \c ODD.)
83
      A = -1,
84
      UNMATCHED = -2  ///< = -2.
76 85
    };
77 86

	
78 87
    typedef typename Graph::template NodeMap<Status> StatusMap;
79 88

	
80 89
  private:
81 90

	
82 91
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
83 92

	
84 93
    typedef UnionFindEnum<IntNodeMap> BlossomSet;
85 94
    typedef ExtendFindEnum<IntNodeMap> TreeSet;
86 95
    typedef RangeMap<Node> NodeIntMap;
87 96
    typedef MatchingMap EarMap;
88 97
    typedef std::vector<Node> NodeQueue;
89 98

	
90 99
    const Graph& _graph;
91 100
    MatchingMap* _matching;
92 101
    StatusMap* _status;
93 102

	
94 103
    EarMap* _ear;
95 104

	
96 105
    IntNodeMap* _blossom_set_index;
97 106
    BlossomSet* _blossom_set;
98 107
    NodeIntMap* _blossom_rep;
99 108

	
100 109
    IntNodeMap* _tree_set_index;
101 110
    TreeSet* _tree_set;
102 111

	
103 112
    NodeQueue _node_queue;
104 113
    int _process, _postpone, _last;
105 114

	
106 115
    int _node_num;
107 116

	
108 117
  private:
109 118

	
110 119
    void createStructures() {
111 120
      _node_num = countNodes(_graph);
112 121
      if (!_matching) {
113 122
        _matching = new MatchingMap(_graph);
114 123
      }
115 124
      if (!_status) {
116 125
        _status = new StatusMap(_graph);
117 126
      }
118 127
      if (!_ear) {
119 128
        _ear = new EarMap(_graph);
120 129
      }
121 130
      if (!_blossom_set) {
122 131
        _blossom_set_index = new IntNodeMap(_graph);
123 132
        _blossom_set = new BlossomSet(*_blossom_set_index);
124 133
      }
125 134
      if (!_blossom_rep) {
126 135
        _blossom_rep = new NodeIntMap(_node_num);
127 136
      }
128 137
      if (!_tree_set) {
129 138
        _tree_set_index = new IntNodeMap(_graph);
130 139
        _tree_set = new TreeSet(*_tree_set_index);
131 140
      }
132 141
      _node_queue.resize(_node_num);
133 142
    }
134 143

	
135 144
    void destroyStructures() {
136 145
      if (_matching) {
137 146
        delete _matching;
138 147
      }
139 148
      if (_status) {
140 149
        delete _status;
141 150
      }
142 151
      if (_ear) {
143 152
        delete _ear;
144 153
      }
145 154
      if (_blossom_set) {
146 155
        delete _blossom_set;
147 156
        delete _blossom_set_index;
148 157
      }
149 158
      if (_blossom_rep) {
150 159
        delete _blossom_rep;
151 160
      }
152 161
      if (_tree_set) {
153 162
        delete _tree_set_index;
154 163
        delete _tree_set;
155 164
      }
156 165
    }
157 166

	
158 167
    void processDense(const Node& n) {
159 168
      _process = _postpone = _last = 0;
160 169
      _node_queue[_last++] = n;
161 170

	
162 171
      while (_process != _last) {
163 172
        Node u = _node_queue[_process++];
164 173
        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
165 174
          Node v = _graph.target(a);
166 175
          if ((*_status)[v] == MATCHED) {
167 176
            extendOnArc(a);
168 177
          } else if ((*_status)[v] == UNMATCHED) {
169 178
            augmentOnArc(a);
170 179
            return;
171 180
          }
172 181
        }
173 182
      }
174 183

	
175 184
      while (_postpone != _last) {
176 185
        Node u = _node_queue[_postpone++];
177 186

	
178 187
        for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
179 188
          Node v = _graph.target(a);
180 189

	
181 190
          if ((*_status)[v] == EVEN) {
182 191
            if (_blossom_set->find(u) != _blossom_set->find(v)) {
183 192
              shrinkOnEdge(a);
184 193
            }
185 194
          }
186 195

	
187 196
          while (_process != _last) {
188 197
            Node w = _node_queue[_process++];
189 198
            for (OutArcIt b(_graph, w); b != INVALID; ++b) {
190 199
              Node x = _graph.target(b);
191 200
              if ((*_status)[x] == MATCHED) {
192 201
                extendOnArc(b);
193 202
              } else if ((*_status)[x] == UNMATCHED) {
194 203
                augmentOnArc(b);
195 204
                return;
196 205
              }
197 206
            }
198 207
          }
199 208
        }
200 209
      }
201 210
    }
202 211

	
203 212
    void processSparse(const Node& n) {
204 213
      _process = _last = 0;
205 214
      _node_queue[_last++] = n;
206 215
      while (_process != _last) {
207 216
        Node u = _node_queue[_process++];
208 217
        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
209 218
          Node v = _graph.target(a);
210 219

	
211 220
          if ((*_status)[v] == EVEN) {
212 221
            if (_blossom_set->find(u) != _blossom_set->find(v)) {
213 222
              shrinkOnEdge(a);
214 223
            }
215 224
          } else if ((*_status)[v] == MATCHED) {
216 225
            extendOnArc(a);
217 226
          } else if ((*_status)[v] == UNMATCHED) {
218 227
            augmentOnArc(a);
219 228
            return;
220 229
          }
221 230
        }
222 231
      }
223 232
    }
224 233

	
225 234
    void shrinkOnEdge(const Edge& e) {
226 235
      Node nca = INVALID;
227 236

	
228 237
      {
229 238
        std::set<Node> left_set, right_set;
230 239

	
231 240
        Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
232 241
        left_set.insert(left);
233 242

	
234 243
        Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
235 244
        right_set.insert(right);
236 245

	
237 246
        while (true) {
238 247
          if ((*_matching)[left] == INVALID) break;
239 248
          left = _graph.target((*_matching)[left]);
240 249
          left = (*_blossom_rep)[_blossom_set->
241 250
                                 find(_graph.target((*_ear)[left]))];
242 251
          if (right_set.find(left) != right_set.end()) {
243 252
            nca = left;
244 253
            break;
245 254
          }
246 255
          left_set.insert(left);
247 256

	
248 257
          if ((*_matching)[right] == INVALID) break;
249 258
          right = _graph.target((*_matching)[right]);
250 259
          right = (*_blossom_rep)[_blossom_set->
251 260
                                  find(_graph.target((*_ear)[right]))];
252 261
          if (left_set.find(right) != left_set.end()) {
253 262
            nca = right;
254 263
            break;
255 264
          }
256 265
          right_set.insert(right);
257 266
        }
258 267

	
259 268
        if (nca == INVALID) {
260 269
          if ((*_matching)[left] == INVALID) {
261 270
            nca = right;
262 271
            while (left_set.find(nca) == left_set.end()) {
263 272
              nca = _graph.target((*_matching)[nca]);
264 273
              nca =(*_blossom_rep)[_blossom_set->
265 274
                                   find(_graph.target((*_ear)[nca]))];
266 275
            }
267 276
          } else {
268 277
            nca = left;
269 278
            while (right_set.find(nca) == right_set.end()) {
270 279
              nca = _graph.target((*_matching)[nca]);
271 280
              nca = (*_blossom_rep)[_blossom_set->
272 281
                                   find(_graph.target((*_ear)[nca]))];
273 282
            }
274 283
          }
275 284
        }
276 285
      }
277 286

	
278 287
      {
279 288

	
280 289
        Node node = _graph.u(e);
281 290
        Arc arc = _graph.direct(e, true);
282 291
        Node base = (*_blossom_rep)[_blossom_set->find(node)];
283 292

	
284 293
        while (base != nca) {
285 294
          (*_ear)[node] = arc;
286 295

	
287 296
          Node n = node;
288 297
          while (n != base) {
289 298
            n = _graph.target((*_matching)[n]);
290 299
            Arc a = (*_ear)[n];
291 300
            n = _graph.target(a);
292 301
            (*_ear)[n] = _graph.oppositeArc(a);
293 302
          }
294 303
          node = _graph.target((*_matching)[base]);
295 304
          _tree_set->erase(base);
296 305
          _tree_set->erase(node);
297 306
          _blossom_set->insert(node, _blossom_set->find(base));
298 307
          (*_status)[node] = EVEN;
299 308
          _node_queue[_last++] = node;
300 309
          arc = _graph.oppositeArc((*_ear)[node]);
301 310
          node = _graph.target((*_ear)[node]);
302 311
          base = (*_blossom_rep)[_blossom_set->find(node)];
303 312
          _blossom_set->join(_graph.target(arc), base);
304 313
        }
305 314
      }
306 315

	
307 316
      (*_blossom_rep)[_blossom_set->find(nca)] = nca;
308 317

	
309 318
      {
310 319

	
311 320
        Node node = _graph.v(e);
312 321
        Arc arc = _graph.direct(e, false);
313 322
        Node base = (*_blossom_rep)[_blossom_set->find(node)];
314 323

	
315 324
        while (base != nca) {
316 325
          (*_ear)[node] = arc;
317 326

	
318 327
          Node n = node;
319 328
          while (n != base) {
320 329
            n = _graph.target((*_matching)[n]);
321 330
            Arc a = (*_ear)[n];
322 331
            n = _graph.target(a);
323 332
            (*_ear)[n] = _graph.oppositeArc(a);
324 333
          }
325 334
          node = _graph.target((*_matching)[base]);
326 335
          _tree_set->erase(base);
327 336
          _tree_set->erase(node);
328 337
          _blossom_set->insert(node, _blossom_set->find(base));
329 338
          (*_status)[node] = EVEN;
330 339
          _node_queue[_last++] = node;
331 340
          arc = _graph.oppositeArc((*_ear)[node]);
332 341
          node = _graph.target((*_ear)[node]);
333 342
          base = (*_blossom_rep)[_blossom_set->find(node)];
334 343
          _blossom_set->join(_graph.target(arc), base);
335 344
        }
336 345
      }
337 346

	
338 347
      (*_blossom_rep)[_blossom_set->find(nca)] = nca;
339 348
    }
340 349

	
341

	
342

	
343 350
    void extendOnArc(const Arc& a) {
344 351
      Node base = _graph.source(a);
345 352
      Node odd = _graph.target(a);
346 353

	
347 354
      (*_ear)[odd] = _graph.oppositeArc(a);
348 355
      Node even = _graph.target((*_matching)[odd]);
349 356
      (*_blossom_rep)[_blossom_set->insert(even)] = even;
350 357
      (*_status)[odd] = ODD;
351 358
      (*_status)[even] = EVEN;
352 359
      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
353 360
      _tree_set->insert(odd, tree);
354 361
      _tree_set->insert(even, tree);
355 362
      _node_queue[_last++] = even;
356 363

	
357 364
    }
358 365

	
359 366
    void augmentOnArc(const Arc& a) {
360 367
      Node even = _graph.source(a);
361 368
      Node odd = _graph.target(a);
362 369

	
363 370
      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
364 371

	
365 372
      (*_matching)[odd] = _graph.oppositeArc(a);
366 373
      (*_status)[odd] = MATCHED;
367 374

	
368 375
      Arc arc = (*_matching)[even];
369 376
      (*_matching)[even] = a;
370 377

	
371 378
      while (arc != INVALID) {
372 379
        odd = _graph.target(arc);
373 380
        arc = (*_ear)[odd];
374 381
        even = _graph.target(arc);
375 382
        (*_matching)[odd] = arc;
376 383
        arc = (*_matching)[even];
377 384
        (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]);
378 385
      }
379 386

	
380 387
      for (typename TreeSet::ItemIt it(*_tree_set, tree);
381 388
           it != INVALID; ++it) {
382 389
        if ((*_status)[it] == ODD) {
383 390
          (*_status)[it] = MATCHED;
384 391
        } else {
385 392
          int blossom = _blossom_set->find(it);
386 393
          for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
387 394
               jt != INVALID; ++jt) {
388 395
            (*_status)[jt] = MATCHED;
389 396
          }
390 397
          _blossom_set->eraseClass(blossom);
391 398
        }
392 399
      }
393 400
      _tree_set->eraseClass(tree);
394 401

	
395 402
    }
396 403

	
397 404
  public:
398 405

	
399 406
    /// \brief Constructor
400 407
    ///
401 408
    /// Constructor.
402 409
    MaxMatching(const Graph& graph)
403 410
      : _graph(graph), _matching(0), _status(0), _ear(0),
404 411
        _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
405 412
        _tree_set_index(0), _tree_set(0) {}
406 413

	
407 414
    ~MaxMatching() {
408 415
      destroyStructures();
409 416
    }
410 417

	
411
    /// \name Execution control
418
    /// \name Execution Control
412 419
    /// The simplest way to execute the algorithm is to use the
413
    /// \c run() member function.
414
    /// \n
415

	
416
    /// If you need better control on the execution, you must call
417
    /// \ref init(), \ref greedyInit() or \ref matchingInit()
418
    /// functions first, then you can start the algorithm with the \ref
419
    /// startSparse() or startDense() functions.
420
    /// \c run() member function.\n
421
    /// If you need better control on the execution, you have to call
422
    /// one of the functions \ref init(), \ref greedyInit() or
423
    /// \ref matchingInit() first, then you can start the algorithm with
424
    /// \ref startSparse() or \ref startDense().
420 425

	
421 426
    ///@{
422 427

	
423
    /// \brief Sets the actual matching to the empty matching.
428
    /// \brief Set the initial matching to the empty matching.
424 429
    ///
425
    /// Sets the actual matching to the empty matching.
426
    ///
430
    /// This function sets the initial matching to the empty matching.
427 431
    void init() {
428 432
      createStructures();
429 433
      for(NodeIt n(_graph); n != INVALID; ++n) {
430 434
        (*_matching)[n] = INVALID;
431 435
        (*_status)[n] = UNMATCHED;
432 436
      }
433 437
    }
434 438

	
435
    ///\brief Finds an initial matching in a greedy way
439
    /// \brief Find an initial matching in a greedy way.
436 440
    ///
437
    ///It finds an initial matching in a greedy way.
441
    /// This function finds an initial matching in a greedy way.
438 442
    void greedyInit() {
439 443
      createStructures();
440 444
      for (NodeIt n(_graph); n != INVALID; ++n) {
441 445
        (*_matching)[n] = INVALID;
442 446
        (*_status)[n] = UNMATCHED;
443 447
      }
444 448
      for (NodeIt n(_graph); n != INVALID; ++n) {
445 449
        if ((*_matching)[n] == INVALID) {
446 450
          for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
447 451
            Node v = _graph.target(a);
448 452
            if ((*_matching)[v] == INVALID && v != n) {
449 453
              (*_matching)[n] = a;
450 454
              (*_status)[n] = MATCHED;
451 455
              (*_matching)[v] = _graph.oppositeArc(a);
452 456
              (*_status)[v] = MATCHED;
453 457
              break;
454 458
            }
455 459
          }
456 460
        }
457 461
      }
458 462
    }
459 463

	
460 464

	
461
    /// \brief Initialize the matching from a map containing.
465
    /// \brief Initialize the matching from a map.
462 466
    ///
463
    /// Initialize the matching from a \c bool valued \c Edge map. This
464
    /// map must have the property that there are no two incident edges
465
    /// with true value, ie. it contains a matching.
467
    /// This function initializes the matching from a \c bool valued edge
468
    /// map. This map should have the property that there are no two incident
469
    /// edges with \c true value, i.e. it really contains a matching.
466 470
    /// \return \c true if the map contains a matching.
467 471
    template <typename MatchingMap>
468 472
    bool matchingInit(const MatchingMap& matching) {
469 473
      createStructures();
470 474

	
471 475
      for (NodeIt n(_graph); n != INVALID; ++n) {
472 476
        (*_matching)[n] = INVALID;
473 477
        (*_status)[n] = UNMATCHED;
474 478
      }
475 479
      for(EdgeIt e(_graph); e!=INVALID; ++e) {
476 480
        if (matching[e]) {
477 481

	
478 482
          Node u = _graph.u(e);
479 483
          if ((*_matching)[u] != INVALID) return false;
480 484
          (*_matching)[u] = _graph.direct(e, true);
481 485
          (*_status)[u] = MATCHED;
482 486

	
483 487
          Node v = _graph.v(e);
484 488
          if ((*_matching)[v] != INVALID) return false;
485 489
          (*_matching)[v] = _graph.direct(e, false);
486 490
          (*_status)[v] = MATCHED;
487 491
        }
488 492
      }
489 493
      return true;
490 494
    }
491 495

	
492
    /// \brief Starts Edmonds' algorithm
496
    /// \brief Start Edmonds' algorithm
493 497
    ///
494
    /// If runs the original Edmonds' algorithm.
498
    /// This function runs the original Edmonds' algorithm.
499
    ///
500
    /// \pre \ref Init(), \ref greedyInit() or \ref matchingInit() must be
501
    /// called before using this function.
495 502
    void startSparse() {
496 503
      for(NodeIt n(_graph); n != INVALID; ++n) {
497 504
        if ((*_status)[n] == UNMATCHED) {
498 505
          (*_blossom_rep)[_blossom_set->insert(n)] = n;
499 506
          _tree_set->insert(n);
500 507
          (*_status)[n] = EVEN;
501 508
          processSparse(n);
502 509
        }
503 510
      }
504 511
    }
505 512

	
506
    /// \brief Starts Edmonds' algorithm.
513
    /// \brief Start Edmonds' algorithm with a heuristic improvement 
514
    /// for dense graphs
507 515
    ///
508
    /// It runs Edmonds' algorithm with a heuristic of postponing
516
    /// This function runs Edmonds' algorithm with a heuristic of postponing
509 517
    /// shrinks, therefore resulting in a faster algorithm for dense graphs.
518
    ///
519
    /// \pre \ref Init(), \ref greedyInit() or \ref matchingInit() must be
520
    /// called before using this function.
510 521
    void startDense() {
511 522
      for(NodeIt n(_graph); n != INVALID; ++n) {
512 523
        if ((*_status)[n] == UNMATCHED) {
513 524
          (*_blossom_rep)[_blossom_set->insert(n)] = n;
514 525
          _tree_set->insert(n);
515 526
          (*_status)[n] = EVEN;
516 527
          processDense(n);
517 528
        }
518 529
      }
519 530
    }
520 531

	
521 532

	
522
    /// \brief Runs Edmonds' algorithm
533
    /// \brief Run Edmonds' algorithm
523 534
    ///
524
    /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
525
    /// or Edmonds' algorithm with a heuristic of
526
    /// postponing shrinks for dense graphs.
535
    /// This function runs Edmonds' algorithm. An additional heuristic of 
536
    /// postponing shrinks is used for relatively dense graphs 
537
    /// (for which <tt>m>=2*n</tt> holds).
527 538
    void run() {
528 539
      if (countEdges(_graph) < 2 * countNodes(_graph)) {
529 540
        greedyInit();
530 541
        startSparse();
531 542
      } else {
532 543
        init();
533 544
        startDense();
534 545
      }
535 546
    }
536 547

	
537 548
    /// @}
538 549

	
539
    /// \name Primal solution
540
    /// Functions to get the primal solution, ie. the matching.
550
    /// \name Primal Solution
551
    /// Functions to get the primal solution, i.e. the maximum matching.
541 552

	
542 553
    /// @{
543 554

	
544
    ///\brief Returns the size of the current matching.
555
    /// \brief Return the size (cardinality) of the matching.
545 556
    ///
546
    ///Returns the size of the current matching. After \ref
547
    ///run() it returns the size of the maximum matching in the graph.
557
    /// This function returns the size (cardinality) of the current matching. 
558
    /// After run() it returns the size of the maximum matching in the graph.
548 559
    int matchingSize() const {
549 560
      int size = 0;
550 561
      for (NodeIt n(_graph); n != INVALID; ++n) {
551 562
        if ((*_matching)[n] != INVALID) {
552 563
          ++size;
553 564
        }
554 565
      }
555 566
      return size / 2;
556 567
    }
557 568

	
558
    /// \brief Returns true when the edge is in the matching.
569
    /// \brief Return \c true if the given edge is in the matching.
559 570
    ///
560
    /// Returns true when the edge is in the matching.
571
    /// This function returns \c true if the given edge is in the current 
572
    /// matching.
561 573
    bool matching(const Edge& edge) const {
562 574
      return edge == (*_matching)[_graph.u(edge)];
563 575
    }
564 576

	
565
    /// \brief Returns the matching edge incident to the given node.
577
    /// \brief Return the matching arc (or edge) incident to the given node.
566 578
    ///
567
    /// Returns the matching edge of a \c node in the actual matching or
568
    /// INVALID if the \c node is not covered by the actual matching.
579
    /// This function returns the matching arc (or edge) incident to the
580
    /// given node in the current matching or \c INVALID if the node is 
581
    /// not covered by the matching.
569 582
    Arc matching(const Node& n) const {
570 583
      return (*_matching)[n];
571 584
    }
572 585

	
573
    ///\brief Returns the mate of a node in the actual matching.
586
    /// \brief Return the mate of the given node.
574 587
    ///
575
    ///Returns the mate of a \c node in the actual matching or
576
    ///INVALID if the \c node is not covered by the actual matching.
588
    /// This function returns the mate of the given node in the current 
589
    /// matching or \c INVALID if the node is not covered by the matching.
577 590
    Node mate(const Node& n) const {
578 591
      return (*_matching)[n] != INVALID ?
579 592
        _graph.target((*_matching)[n]) : INVALID;
580 593
    }
581 594

	
582 595
    /// @}
583 596

	
584
    /// \name Dual solution
585
    /// Functions to get the dual solution, ie. the decomposition.
597
    /// \name Dual Solution
598
    /// Functions to get the dual solution, i.e. the Gallai-Edmonds 
599
    /// decomposition.
586 600

	
587 601
    /// @{
588 602

	
589
    /// \brief Returns the class of the node in the Edmonds-Gallai
603
    /// \brief Return the status of the given node in the Edmonds-Gallai
590 604
    /// decomposition.
591 605
    ///
592
    /// Returns the class of the node in the Edmonds-Gallai
593
    /// decomposition.
606
    /// This function returns the \ref Status "status" of the given node
607
    /// in the Edmonds-Gallai decomposition.
594 608
    Status decomposition(const Node& n) const {
595 609
      return (*_status)[n];
596 610
    }
597 611

	
598
    /// \brief Returns true when the node is in the barrier.
612
    /// \brief Return \c true if the given node is in the barrier.
599 613
    ///
600
    /// Returns true when the node is in the barrier.
614
    /// This function returns \c true if the given node is in the barrier.
601 615
    bool barrier(const Node& n) const {
602 616
      return (*_status)[n] == ODD;
603 617
    }
604 618

	
605 619
    /// @}
606 620

	
607 621
  };
608 622

	
609 623
  /// \ingroup matching
610 624
  ///
611 625
  /// \brief Weighted matching in general graphs
612 626
  ///
613 627
  /// This class provides an efficient implementation of Edmond's
614 628
  /// maximum weighted matching algorithm. The implementation is based
615 629
  /// on extensive use of priority queues and provides
616 630
  /// \f$O(nm\log n)\f$ time complexity.
617 631
  ///
618
  /// The maximum weighted matching problem is to find undirected
619
  /// edges in the graph with maximum overall weight and no two of
620
  /// them shares their ends. The problem can be formulated with the
621
  /// following linear program.
632
  /// The maximum weighted matching problem is to find a subset of the 
633
  /// edges in an undirected graph with maximum overall weight for which 
634
  /// each node has at most one incident edge.
635
  /// It can be formulated with the following linear program.
622 636
  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
623 637
  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
624 638
      \quad \forall B\in\mathcal{O}\f] */
625 639
  /// \f[x_e \ge 0\quad \forall e\in E\f]
626 640
  /// \f[\max \sum_{e\in E}x_ew_e\f]
627 641
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
628 642
  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
629 643
  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
630 644
  /// subsets of the nodes.
631 645
  ///
632 646
  /// The algorithm calculates an optimal matching and a proof of the
633 647
  /// optimality. The solution of the dual problem can be used to check
634 648
  /// the result of the algorithm. The dual linear problem is the
649
  /// following.
635 650
  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
636 651
      z_B \ge w_{uv} \quad \forall uv\in E\f] */
637 652
  /// \f[y_u \ge 0 \quad \forall u \in V\f]
638 653
  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
639 654
  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
640 655
      \frac{\vert B \vert - 1}{2}z_B\f] */
641 656
  ///
642
  /// The algorithm can be executed with \c run() or the \c init() and
643
  /// then the \c start() member functions. After it the matching can
644
  /// be asked with \c matching() or mate() functions. The dual
645
  /// solution can be get with \c nodeValue(), \c blossomNum() and \c
646
  /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
647
  /// "BlossomIt" nested class, which is able to iterate on the nodes
648
  /// of a blossom. If the value type is integral then the dual
649
  /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
657
  /// The algorithm can be executed with the run() function. 
658
  /// After it the matching (the primal solution) and the dual solution
659
  /// can be obtained using the query functions and the 
660
  /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, 
661
  /// which is able to iterate on the nodes of a blossom. 
662
  /// If the value type is integer, then the dual solution is multiplied
663
  /// by \ref MaxWeightedMatching::dualScale "4".
664
  ///
665
  /// \tparam GR The graph type the algorithm runs on.
666
  /// \tparam WM The type edge weight map. The default type is 
667
  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
668
#ifdef DOXYGEN
669
  template <typename GR, typename WM>
670
#else
650 671
  template <typename GR,
651 672
            typename WM = typename GR::template EdgeMap<int> >
673
#endif
652 674
  class MaxWeightedMatching {
653 675
  public:
654 676

	
655
    ///\e
677
    /// The graph type of the algorithm
656 678
    typedef GR Graph;
657
    ///\e
679
    /// The type of the edge weight map
658 680
    typedef WM WeightMap;
659
    ///\e
681
    /// The value type of the edge weights
660 682
    typedef typename WeightMap::Value Value;
661 683

	
684
    typedef typename Graph::template NodeMap<typename Graph::Arc>
685
    MatchingMap;
686

	
662 687
    /// \brief Scaling factor for dual solution
663 688
    ///
664
    /// Scaling factor for dual solution, it is equal to 4 or 1
689
    /// Scaling factor for dual solution. It is equal to 4 or 1
665 690
    /// according to the value type.
666 691
    static const int dualScale =
667 692
      std::numeric_limits<Value>::is_integer ? 4 : 1;
668 693

	
669
    typedef typename Graph::template NodeMap<typename Graph::Arc>
670
    MatchingMap;
671

	
672 694
  private:
673 695

	
674 696
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
675 697

	
676 698
    typedef typename Graph::template NodeMap<Value> NodePotential;
677 699
    typedef std::vector<Node> BlossomNodeList;
678 700

	
679 701
    struct BlossomVariable {
680 702
      int begin, end;
681 703
      Value value;
682 704

	
683 705
      BlossomVariable(int _begin, int _end, Value _value)
684 706
        : begin(_begin), end(_end), value(_value) {}
685 707

	
686 708
    };
687 709

	
688 710
    typedef std::vector<BlossomVariable> BlossomPotential;
689 711

	
690 712
    const Graph& _graph;
691 713
    const WeightMap& _weight;
692 714

	
693 715
    MatchingMap* _matching;
694 716

	
695 717
    NodePotential* _node_potential;
696 718

	
697 719
    BlossomPotential _blossom_potential;
698 720
    BlossomNodeList _blossom_node_list;
699 721

	
700 722
    int _node_num;
701 723
    int _blossom_num;
702 724

	
703 725
    typedef RangeMap<int> IntIntMap;
704 726

	
705 727
    enum Status {
706 728
      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
707 729
    };
708 730

	
709 731
    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
710 732
    struct BlossomData {
711 733
      int tree;
712 734
      Status status;
713 735
      Arc pred, next;
714 736
      Value pot, offset;
715 737
      Node base;
716 738
    };
717 739

	
718 740
    IntNodeMap *_blossom_index;
719 741
    BlossomSet *_blossom_set;
720 742
    RangeMap<BlossomData>* _blossom_data;
721 743

	
722 744
    IntNodeMap *_node_index;
723 745
    IntArcMap *_node_heap_index;
724 746

	
725 747
    struct NodeData {
726 748

	
727 749
      NodeData(IntArcMap& node_heap_index)
728 750
        : heap(node_heap_index) {}
729 751

	
730 752
      int blossom;
731 753
      Value pot;
732 754
      BinHeap<Value, IntArcMap> heap;
733 755
      std::map<int, Arc> heap_index;
734 756

	
735 757
      int tree;
736 758
    };
737 759

	
738 760
    RangeMap<NodeData>* _node_data;
739 761

	
740 762
    typedef ExtendFindEnum<IntIntMap> TreeSet;
741 763

	
742 764
    IntIntMap *_tree_set_index;
743 765
    TreeSet *_tree_set;
744 766

	
745 767
    IntNodeMap *_delta1_index;
746 768
    BinHeap<Value, IntNodeMap> *_delta1;
747 769

	
748 770
    IntIntMap *_delta2_index;
749 771
    BinHeap<Value, IntIntMap> *_delta2;
750 772

	
751 773
    IntEdgeMap *_delta3_index;
752 774
    BinHeap<Value, IntEdgeMap> *_delta3;
753 775

	
754 776
    IntIntMap *_delta4_index;
755 777
    BinHeap<Value, IntIntMap> *_delta4;
756 778

	
757 779
    Value _delta_sum;
758 780

	
759 781
    void createStructures() {
760 782
      _node_num = countNodes(_graph);
761 783
      _blossom_num = _node_num * 3 / 2;
762 784

	
763 785
      if (!_matching) {
764 786
        _matching = new MatchingMap(_graph);
765 787
      }
766 788
      if (!_node_potential) {
767 789
        _node_potential = new NodePotential(_graph);
768 790
      }
769 791
      if (!_blossom_set) {
770 792
        _blossom_index = new IntNodeMap(_graph);
771 793
        _blossom_set = new BlossomSet(*_blossom_index);
772 794
        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
773 795
      }
774 796

	
775 797
      if (!_node_index) {
776 798
        _node_index = new IntNodeMap(_graph);
777 799
        _node_heap_index = new IntArcMap(_graph);
778 800
        _node_data = new RangeMap<NodeData>(_node_num,
779 801
                                              NodeData(*_node_heap_index));
780 802
      }
781 803

	
782 804
      if (!_tree_set) {
783 805
        _tree_set_index = new IntIntMap(_blossom_num);
784 806
        _tree_set = new TreeSet(*_tree_set_index);
785 807
      }
786 808
      if (!_delta1) {
787 809
        _delta1_index = new IntNodeMap(_graph);
788 810
        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
789 811
      }
790 812
      if (!_delta2) {
791 813
        _delta2_index = new IntIntMap(_blossom_num);
792 814
        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
793 815
      }
794 816
      if (!_delta3) {
795 817
        _delta3_index = new IntEdgeMap(_graph);
796 818
        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
797 819
      }
798 820
      if (!_delta4) {
799 821
        _delta4_index = new IntIntMap(_blossom_num);
800 822
        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
801 823
      }
802 824
    }
803 825

	
804 826
    void destroyStructures() {
805 827
      _node_num = countNodes(_graph);
806 828
      _blossom_num = _node_num * 3 / 2;
807 829

	
808 830
      if (_matching) {
809 831
        delete _matching;
810 832
      }
811 833
      if (_node_potential) {
812 834
        delete _node_potential;
813 835
      }
814 836
      if (_blossom_set) {
815 837
        delete _blossom_index;
816 838
        delete _blossom_set;
817 839
        delete _blossom_data;
818 840
      }
819 841

	
820 842
      if (_node_index) {
821 843
        delete _node_index;
822 844
        delete _node_heap_index;
823 845
        delete _node_data;
824 846
      }
825 847

	
826 848
      if (_tree_set) {
827 849
        delete _tree_set_index;
828 850
        delete _tree_set;
829 851
      }
830 852
      if (_delta1) {
831 853
        delete _delta1_index;
832 854
        delete _delta1;
833 855
      }
834 856
      if (_delta2) {
835 857
        delete _delta2_index;
836 858
        delete _delta2;
837 859
      }
838 860
      if (_delta3) {
839 861
        delete _delta3_index;
840 862
        delete _delta3;
841 863
      }
842 864
      if (_delta4) {
843 865
        delete _delta4_index;
844 866
        delete _delta4;
845 867
      }
846 868
    }
847 869

	
848 870
    void matchedToEven(int blossom, int tree) {
849 871
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
850 872
        _delta2->erase(blossom);
851 873
      }
852 874

	
853 875
      if (!_blossom_set->trivial(blossom)) {
854 876
        (*_blossom_data)[blossom].pot -=
855 877
          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
856 878
      }
857 879

	
858 880
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
859 881
           n != INVALID; ++n) {
860 882

	
861 883
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
862 884
        int ni = (*_node_index)[n];
863 885

	
... ...
@@ -1442,753 +1464,798 @@
1442 1464
      oddToMatched(blossom);
1443 1465
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1444 1466
        _delta2->erase(blossom);
1445 1467
      }
1446 1468

	
1447 1469
      std::vector<int> subblossoms;
1448 1470
      _blossom_set->split(blossom, std::back_inserter(subblossoms));
1449 1471

	
1450 1472
      Value offset = (*_blossom_data)[blossom].offset;
1451 1473
      int b = _blossom_set->find(_graph.source(pred));
1452 1474
      int d = _blossom_set->find(_graph.source(next));
1453 1475

	
1454 1476
      int ib = -1, id = -1;
1455 1477
      for (int i = 0; i < int(subblossoms.size()); ++i) {
1456 1478
        if (subblossoms[i] == b) ib = i;
1457 1479
        if (subblossoms[i] == d) id = i;
1458 1480

	
1459 1481
        (*_blossom_data)[subblossoms[i]].offset = offset;
1460 1482
        if (!_blossom_set->trivial(subblossoms[i])) {
1461 1483
          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1462 1484
        }
1463 1485
        if (_blossom_set->classPrio(subblossoms[i]) !=
1464 1486
            std::numeric_limits<Value>::max()) {
1465 1487
          _delta2->push(subblossoms[i],
1466 1488
                        _blossom_set->classPrio(subblossoms[i]) -
1467 1489
                        (*_blossom_data)[subblossoms[i]].offset);
1468 1490
        }
1469 1491
      }
1470 1492

	
1471 1493
      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1472 1494
        for (int i = (id + 1) % subblossoms.size();
1473 1495
             i != ib; i = (i + 2) % subblossoms.size()) {
1474 1496
          int sb = subblossoms[i];
1475 1497
          int tb = subblossoms[(i + 1) % subblossoms.size()];
1476 1498
          (*_blossom_data)[sb].next =
1477 1499
            _graph.oppositeArc((*_blossom_data)[tb].next);
1478 1500
        }
1479 1501

	
1480 1502
        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1481 1503
          int sb = subblossoms[i];
1482 1504
          int tb = subblossoms[(i + 1) % subblossoms.size()];
1483 1505
          int ub = subblossoms[(i + 2) % subblossoms.size()];
1484 1506

	
1485 1507
          (*_blossom_data)[sb].status = ODD;
1486 1508
          matchedToOdd(sb);
1487 1509
          _tree_set->insert(sb, tree);
1488 1510
          (*_blossom_data)[sb].pred = pred;
1489 1511
          (*_blossom_data)[sb].next =
1490 1512
                           _graph.oppositeArc((*_blossom_data)[tb].next);
1491 1513

	
1492 1514
          pred = (*_blossom_data)[ub].next;
1493 1515

	
1494 1516
          (*_blossom_data)[tb].status = EVEN;
1495 1517
          matchedToEven(tb, tree);
1496 1518
          _tree_set->insert(tb, tree);
1497 1519
          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1498 1520
        }
1499 1521

	
1500 1522
        (*_blossom_data)[subblossoms[id]].status = ODD;
1501 1523
        matchedToOdd(subblossoms[id]);
1502 1524
        _tree_set->insert(subblossoms[id], tree);
1503 1525
        (*_blossom_data)[subblossoms[id]].next = next;
1504 1526
        (*_blossom_data)[subblossoms[id]].pred = pred;
1505 1527

	
1506 1528
      } else {
1507 1529

	
1508 1530
        for (int i = (ib + 1) % subblossoms.size();
1509 1531
             i != id; i = (i + 2) % subblossoms.size()) {
1510 1532
          int sb = subblossoms[i];
1511 1533
          int tb = subblossoms[(i + 1) % subblossoms.size()];
1512 1534
          (*_blossom_data)[sb].next =
1513 1535
            _graph.oppositeArc((*_blossom_data)[tb].next);
1514 1536
        }
1515 1537

	
1516 1538
        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1517 1539
          int sb = subblossoms[i];
1518 1540
          int tb = subblossoms[(i + 1) % subblossoms.size()];
1519 1541
          int ub = subblossoms[(i + 2) % subblossoms.size()];
1520 1542

	
1521 1543
          (*_blossom_data)[sb].status = ODD;
1522 1544
          matchedToOdd(sb);
1523 1545
          _tree_set->insert(sb, tree);
1524 1546
          (*_blossom_data)[sb].next = next;
1525 1547
          (*_blossom_data)[sb].pred =
1526 1548
            _graph.oppositeArc((*_blossom_data)[tb].next);
1527 1549

	
1528 1550
          (*_blossom_data)[tb].status = EVEN;
1529 1551
          matchedToEven(tb, tree);
1530 1552
          _tree_set->insert(tb, tree);
1531 1553
          (*_blossom_data)[tb].pred =
1532 1554
            (*_blossom_data)[tb].next =
1533 1555
            _graph.oppositeArc((*_blossom_data)[ub].next);
1534 1556
          next = (*_blossom_data)[ub].next;
1535 1557
        }
1536 1558

	
1537 1559
        (*_blossom_data)[subblossoms[ib]].status = ODD;
1538 1560
        matchedToOdd(subblossoms[ib]);
1539 1561
        _tree_set->insert(subblossoms[ib], tree);
1540 1562
        (*_blossom_data)[subblossoms[ib]].next = next;
1541 1563
        (*_blossom_data)[subblossoms[ib]].pred = pred;
1542 1564
      }
1543 1565
      _tree_set->erase(blossom);
1544 1566
    }
1545 1567

	
1546 1568
    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1547 1569
      if (_blossom_set->trivial(blossom)) {
1548 1570
        int bi = (*_node_index)[base];
1549 1571
        Value pot = (*_node_data)[bi].pot;
1550 1572

	
1551 1573
        (*_matching)[base] = matching;
1552 1574
        _blossom_node_list.push_back(base);
1553 1575
        (*_node_potential)[base] = pot;
1554 1576
      } else {
1555 1577

	
1556 1578
        Value pot = (*_blossom_data)[blossom].pot;
1557 1579
        int bn = _blossom_node_list.size();
1558 1580

	
1559 1581
        std::vector<int> subblossoms;
1560 1582
        _blossom_set->split(blossom, std::back_inserter(subblossoms));
1561 1583
        int b = _blossom_set->find(base);
1562 1584
        int ib = -1;
1563 1585
        for (int i = 0; i < int(subblossoms.size()); ++i) {
1564 1586
          if (subblossoms[i] == b) { ib = i; break; }
1565 1587
        }
1566 1588

	
1567 1589
        for (int i = 1; i < int(subblossoms.size()); i += 2) {
1568 1590
          int sb = subblossoms[(ib + i) % subblossoms.size()];
1569 1591
          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1570 1592

	
1571 1593
          Arc m = (*_blossom_data)[tb].next;
1572 1594
          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1573 1595
          extractBlossom(tb, _graph.source(m), m);
1574 1596
        }
1575 1597
        extractBlossom(subblossoms[ib], base, matching);
1576 1598

	
1577 1599
        int en = _blossom_node_list.size();
1578 1600

	
1579 1601
        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1580 1602
      }
1581 1603
    }
1582 1604

	
1583 1605
    void extractMatching() {
1584 1606
      std::vector<int> blossoms;
1585 1607
      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1586 1608
        blossoms.push_back(c);
1587 1609
      }
1588 1610

	
1589 1611
      for (int i = 0; i < int(blossoms.size()); ++i) {
1590 1612
        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
1591 1613

	
1592 1614
          Value offset = (*_blossom_data)[blossoms[i]].offset;
1593 1615
          (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1594 1616
          for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1595 1617
               n != INVALID; ++n) {
1596 1618
            (*_node_data)[(*_node_index)[n]].pot -= offset;
1597 1619
          }
1598 1620

	
1599 1621
          Arc matching = (*_blossom_data)[blossoms[i]].next;
1600 1622
          Node base = _graph.source(matching);
1601 1623
          extractBlossom(blossoms[i], base, matching);
1602 1624
        } else {
1603 1625
          Node base = (*_blossom_data)[blossoms[i]].base;
1604 1626
          extractBlossom(blossoms[i], base, INVALID);
1605 1627
        }
1606 1628
      }
1607 1629
    }
1608 1630

	
1609 1631
  public:
1610 1632

	
1611 1633
    /// \brief Constructor
1612 1634
    ///
1613 1635
    /// Constructor.
1614 1636
    MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1615 1637
      : _graph(graph), _weight(weight), _matching(0),
1616 1638
        _node_potential(0), _blossom_potential(), _blossom_node_list(),
1617 1639
        _node_num(0), _blossom_num(0),
1618 1640

	
1619 1641
        _blossom_index(0), _blossom_set(0), _blossom_data(0),
1620 1642
        _node_index(0), _node_heap_index(0), _node_data(0),
1621 1643
        _tree_set_index(0), _tree_set(0),
1622 1644

	
1623 1645
        _delta1_index(0), _delta1(0),
1624 1646
        _delta2_index(0), _delta2(0),
1625 1647
        _delta3_index(0), _delta3(0),
1626 1648
        _delta4_index(0), _delta4(0),
1627 1649

	
1628 1650
        _delta_sum() {}
1629 1651

	
1630 1652
    ~MaxWeightedMatching() {
1631 1653
      destroyStructures();
1632 1654
    }
1633 1655

	
1634
    /// \name Execution control
1656
    /// \name Execution Control
1635 1657
    /// The simplest way to execute the algorithm is to use the
1636
    /// \c run() member function.
1658
    /// \ref run() member function.
1637 1659

	
1638 1660
    ///@{
1639 1661

	
1640 1662
    /// \brief Initialize the algorithm
1641 1663
    ///
1642
    /// Initialize the algorithm
1664
    /// This function initializes the algorithm.
1643 1665
    void init() {
1644 1666
      createStructures();
1645 1667

	
1646 1668
      for (ArcIt e(_graph); e != INVALID; ++e) {
1647 1669
        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1648 1670
      }
1649 1671
      for (NodeIt n(_graph); n != INVALID; ++n) {
1650 1672
        (*_delta1_index)[n] = _delta1->PRE_HEAP;
1651 1673
      }
1652 1674
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1653 1675
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1654 1676
      }
1655 1677
      for (int i = 0; i < _blossom_num; ++i) {
1656 1678
        (*_delta2_index)[i] = _delta2->PRE_HEAP;
1657 1679
        (*_delta4_index)[i] = _delta4->PRE_HEAP;
1658 1680
      }
1659 1681

	
1660 1682
      int index = 0;
1661 1683
      for (NodeIt n(_graph); n != INVALID; ++n) {
1662 1684
        Value max = 0;
1663 1685
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1664 1686
          if (_graph.target(e) == n) continue;
1665 1687
          if ((dualScale * _weight[e]) / 2 > max) {
1666 1688
            max = (dualScale * _weight[e]) / 2;
1667 1689
          }
1668 1690
        }
1669 1691
        (*_node_index)[n] = index;
1670 1692
        (*_node_data)[index].pot = max;
1671 1693
        _delta1->push(n, max);
1672 1694
        int blossom =
1673 1695
          _blossom_set->insert(n, std::numeric_limits<Value>::max());
1674 1696

	
1675 1697
        _tree_set->insert(blossom);
1676 1698

	
1677 1699
        (*_blossom_data)[blossom].status = EVEN;
1678 1700
        (*_blossom_data)[blossom].pred = INVALID;
1679 1701
        (*_blossom_data)[blossom].next = INVALID;
1680 1702
        (*_blossom_data)[blossom].pot = 0;
1681 1703
        (*_blossom_data)[blossom].offset = 0;
1682 1704
        ++index;
1683 1705
      }
1684 1706
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1685 1707
        int si = (*_node_index)[_graph.u(e)];
1686 1708
        int ti = (*_node_index)[_graph.v(e)];
1687 1709
        if (_graph.u(e) != _graph.v(e)) {
1688 1710
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1689 1711
                            dualScale * _weight[e]) / 2);
1690 1712
        }
1691 1713
      }
1692 1714
    }
1693 1715

	
1694
    /// \brief Starts the algorithm
1716
    /// \brief Start the algorithm
1695 1717
    ///
1696
    /// Starts the algorithm
1718
    /// This function starts the algorithm.
1719
    ///
1720
    /// \pre \ref init() must be called before using this function.
1697 1721
    void start() {
1698 1722
      enum OpType {
1699 1723
        D1, D2, D3, D4
1700 1724
      };
1701 1725

	
1702 1726
      int unmatched = _node_num;
1703 1727
      while (unmatched > 0) {
1704 1728
        Value d1 = !_delta1->empty() ?
1705 1729
          _delta1->prio() : std::numeric_limits<Value>::max();
1706 1730

	
1707 1731
        Value d2 = !_delta2->empty() ?
1708 1732
          _delta2->prio() : std::numeric_limits<Value>::max();
1709 1733

	
1710 1734
        Value d3 = !_delta3->empty() ?
1711 1735
          _delta3->prio() : std::numeric_limits<Value>::max();
1712 1736

	
1713 1737
        Value d4 = !_delta4->empty() ?
1714 1738
          _delta4->prio() : std::numeric_limits<Value>::max();
1715 1739

	
1716 1740
        _delta_sum = d1; OpType ot = D1;
1717 1741
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1718 1742
        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1719 1743
        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1720 1744

	
1721 1745

	
1722 1746
        switch (ot) {
1723 1747
        case D1:
1724 1748
          {
1725 1749
            Node n = _delta1->top();
1726 1750
            unmatchNode(n);
1727 1751
            --unmatched;
1728 1752
          }
1729 1753
          break;
1730 1754
        case D2:
1731 1755
          {
1732 1756
            int blossom = _delta2->top();
1733 1757
            Node n = _blossom_set->classTop(blossom);
1734 1758
            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
1735 1759
            extendOnArc(e);
1736 1760
          }
1737 1761
          break;
1738 1762
        case D3:
1739 1763
          {
1740 1764
            Edge e = _delta3->top();
1741 1765

	
1742 1766
            int left_blossom = _blossom_set->find(_graph.u(e));
1743 1767
            int right_blossom = _blossom_set->find(_graph.v(e));
1744 1768

	
1745 1769
            if (left_blossom == right_blossom) {
1746 1770
              _delta3->pop();
1747 1771
            } else {
1748 1772
              int left_tree;
1749 1773
              if ((*_blossom_data)[left_blossom].status == EVEN) {
1750 1774
                left_tree = _tree_set->find(left_blossom);
1751 1775
              } else {
1752 1776
                left_tree = -1;
1753 1777
                ++unmatched;
1754 1778
              }
1755 1779
              int right_tree;
1756 1780
              if ((*_blossom_data)[right_blossom].status == EVEN) {
1757 1781
                right_tree = _tree_set->find(right_blossom);
1758 1782
              } else {
1759 1783
                right_tree = -1;
1760 1784
                ++unmatched;
1761 1785
              }
1762 1786

	
1763 1787
              if (left_tree == right_tree) {
1764 1788
                shrinkOnEdge(e, left_tree);
1765 1789
              } else {
1766 1790
                augmentOnEdge(e);
1767 1791
                unmatched -= 2;
1768 1792
              }
1769 1793
            }
1770 1794
          } break;
1771 1795
        case D4:
1772 1796
          splitBlossom(_delta4->top());
1773 1797
          break;
1774 1798
        }
1775 1799
      }
1776 1800
      extractMatching();
1777 1801
    }
1778 1802

	
1779
    /// \brief Runs %MaxWeightedMatching algorithm.
1803
    /// \brief Run the algorithm.
1780 1804
    ///
1781
    /// This method runs the %MaxWeightedMatching algorithm.
1805
    /// This method runs the \c %MaxWeightedMatching algorithm.
1782 1806
    ///
1783 1807
    /// \note mwm.run() is just a shortcut of the following code.
1784 1808
    /// \code
1785 1809
    ///   mwm.init();
1786 1810
    ///   mwm.start();
1787 1811
    /// \endcode
1788 1812
    void run() {
1789 1813
      init();
1790 1814
      start();
1791 1815
    }
1792 1816

	
1793 1817
    /// @}
1794 1818

	
1795
    /// \name Primal solution
1796
    /// Functions to get the primal solution, ie. the matching.
1819
    /// \name Primal Solution
1820
    /// Functions to get the primal solution, i.e. the maximum weighted 
1821
    /// matching.\n
1822
    /// Either \ref run() or \ref start() function should be called before
1823
    /// using them.
1797 1824

	
1798 1825
    /// @{
1799 1826

	
1800
    /// \brief Returns the weight of the matching.
1827
    /// \brief Return the weight of the matching.
1801 1828
    ///
1802
    /// Returns the weight of the matching.
1829
    /// This function returns the weight of the found matching.
1830
    ///
1831
    /// \pre Either run() or start() must be called before using this function.
1803 1832
    Value matchingValue() const {
1804 1833
      Value sum = 0;
1805 1834
      for (NodeIt n(_graph); n != INVALID; ++n) {
1806 1835
        if ((*_matching)[n] != INVALID) {
1807 1836
          sum += _weight[(*_matching)[n]];
1808 1837
        }
1809 1838
      }
1810 1839
      return sum /= 2;
1811 1840
    }
1812 1841

	
1813
    /// \brief Returns the cardinality of the matching.
1842
    /// \brief Return the size (cardinality) of the matching.
1814 1843
    ///
1815
    /// Returns the cardinality of the matching.
1844
    /// This function returns the size (cardinality) of the found matching.
1845
    ///
1846
    /// \pre Either run() or start() must be called before using this function.
1816 1847
    int matchingSize() const {
1817 1848
      int num = 0;
1818 1849
      for (NodeIt n(_graph); n != INVALID; ++n) {
1819 1850
        if ((*_matching)[n] != INVALID) {
1820 1851
          ++num;
1821 1852
        }
1822 1853
      }
1823 1854
      return num /= 2;
1824 1855
    }
1825 1856

	
1826
    /// \brief Returns true when the edge is in the matching.
1857
    /// \brief Return \c true if the given edge is in the matching.
1827 1858
    ///
1828
    /// Returns true when the edge is in the matching.
1859
    /// This function returns \c true if the given edge is in the found 
1860
    /// matching.
1861
    ///
1862
    /// \pre Either run() or start() must be called before using this function.
1829 1863
    bool matching(const Edge& edge) const {
1830 1864
      return edge == (*_matching)[_graph.u(edge)];
1831 1865
    }
1832 1866

	
1833
    /// \brief Returns the incident matching arc.
1867
    /// \brief Return the matching arc (or edge) incident to the given node.
1834 1868
    ///
1835
    /// Returns the incident matching arc from given node. If the
1836
    /// node is not matched then it gives back \c INVALID.
1869
    /// This function returns the matching arc (or edge) incident to the
1870
    /// given node in the found matching or \c INVALID if the node is 
1871
    /// not covered by the matching.
1872
    ///
1873
    /// \pre Either run() or start() must be called before using this function.
1837 1874
    Arc matching(const Node& node) const {
1838 1875
      return (*_matching)[node];
1839 1876
    }
1840 1877

	
1841
    /// \brief Returns the mate of the node.
1878
    /// \brief Return the mate of the given node.
1842 1879
    ///
1843
    /// Returns the adjancent node in a mathcing arc. If the node is
1844
    /// not matched then it gives back \c INVALID.
1880
    /// This function returns the mate of the given node in the found 
1881
    /// matching or \c INVALID if the node is not covered by the matching.
1882
    ///
1883
    /// \pre Either run() or start() must be called before using this function.
1845 1884
    Node mate(const Node& node) const {
1846 1885
      return (*_matching)[node] != INVALID ?
1847 1886
        _graph.target((*_matching)[node]) : INVALID;
1848 1887
    }
1849 1888

	
1850 1889
    /// @}
1851 1890

	
1852
    /// \name Dual solution
1853
    /// Functions to get the dual solution.
1891
    /// \name Dual Solution
1892
    /// Functions to get the dual solution.\n
1893
    /// Either \ref run() or \ref start() function should be called before
1894
    /// using them.
1854 1895

	
1855 1896
    /// @{
1856 1897

	
1857
    /// \brief Returns the value of the dual solution.
1898
    /// \brief Return the value of the dual solution.
1858 1899
    ///
1859
    /// Returns the value of the dual solution. It should be equal to
1860
    /// the primal value scaled by \ref dualScale "dual scale".
1900
    /// This function returns the value of the dual solution. 
1901
    /// It should be equal to the primal value scaled by \ref dualScale 
1902
    /// "dual scale".
1903
    ///
1904
    /// \pre Either run() or start() must be called before using this function.
1861 1905
    Value dualValue() const {
1862 1906
      Value sum = 0;
1863 1907
      for (NodeIt n(_graph); n != INVALID; ++n) {
1864 1908
        sum += nodeValue(n);
1865 1909
      }
1866 1910
      for (int i = 0; i < blossomNum(); ++i) {
1867 1911
        sum += blossomValue(i) * (blossomSize(i) / 2);
1868 1912
      }
1869 1913
      return sum;
1870 1914
    }
1871 1915

	
1872
    /// \brief Returns the value of the node.
1916
    /// \brief Return the dual value (potential) of the given node.
1873 1917
    ///
1874
    /// Returns the the value of the node.
1918
    /// This function returns the dual value (potential) of the given node.
1919
    ///
1920
    /// \pre Either run() or start() must be called before using this function.
1875 1921
    Value nodeValue(const Node& n) const {
1876 1922
      return (*_node_potential)[n];
1877 1923
    }
1878 1924

	
1879
    /// \brief Returns the number of the blossoms in the basis.
1925
    /// \brief Return the number of the blossoms in the basis.
1880 1926
    ///
1881
    /// Returns the number of the blossoms in the basis.
1927
    /// This function returns the number of the blossoms in the basis.
1928
    ///
1929
    /// \pre Either run() or start() must be called before using this function.
1882 1930
    /// \see BlossomIt
1883 1931
    int blossomNum() const {
1884 1932
      return _blossom_potential.size();
1885 1933
    }
1886 1934

	
1887

	
1888
    /// \brief Returns the number of the nodes in the blossom.
1935
    /// \brief Return the number of the nodes in the given blossom.
1889 1936
    ///
1890
    /// Returns the number of the nodes in the blossom.
1937
    /// This function returns the number of the nodes in the given blossom.
1938
    ///
1939
    /// \pre Either run() or start() must be called before using this function.
1940
    /// \see BlossomIt
1891 1941
    int blossomSize(int k) const {
1892 1942
      return _blossom_potential[k].end - _blossom_potential[k].begin;
1893 1943
    }
1894 1944

	
1895
    /// \brief Returns the value of the blossom.
1945
    /// \brief Return the dual value (ptential) of the given blossom.
1896 1946
    ///
1897
    /// Returns the the value of the blossom.
1898
    /// \see BlossomIt
1947
    /// This function returns the dual value (ptential) of the given blossom.
1948
    ///
1949
    /// \pre Either run() or start() must be called before using this function.
1899 1950
    Value blossomValue(int k) const {
1900 1951
      return _blossom_potential[k].value;
1901 1952
    }
1902 1953

	
1903
    /// \brief Iterator for obtaining the nodes of the blossom.
1954
    /// \brief Iterator for obtaining the nodes of a blossom.
1904 1955
    ///
1905
    /// Iterator for obtaining the nodes of the blossom. This class
1906
    /// provides a common lemon style iterator for listing a
1907
    /// subset of the nodes.
1956
    /// This class provides an iterator for obtaining the nodes of the 
1957
    /// given blossom. It lists a subset of the nodes.
1958
    /// Before using this iterator, you must allocate a 
1959
    /// MaxWeightedMatching class and execute it.
1908 1960
    class BlossomIt {
1909 1961
    public:
1910 1962

	
1911 1963
      /// \brief Constructor.
1912 1964
      ///
1913
      /// Constructor to get the nodes of the variable.
1965
      /// Constructor to get the nodes of the given variable.
1966
      ///
1967
      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 
1968
      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 
1969
      /// called before initializing this iterator.
1914 1970
      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
1915 1971
        : _algorithm(&algorithm)
1916 1972
      {
1917 1973
        _index = _algorithm->_blossom_potential[variable].begin;
1918 1974
        _last = _algorithm->_blossom_potential[variable].end;
1919 1975
      }
1920 1976

	
1921
      /// \brief Conversion to node.
1977
      /// \brief Conversion to \c Node.
1922 1978
      ///
1923
      /// Conversion to node.
1979
      /// Conversion to \c Node.
1924 1980
      operator Node() const {
1925 1981
        return _algorithm->_blossom_node_list[_index];
1926 1982
      }
1927 1983

	
1928 1984
      /// \brief Increment operator.
1929 1985
      ///
1930 1986
      /// Increment operator.
1931 1987
      BlossomIt& operator++() {
1932 1988
        ++_index;
1933 1989
        return *this;
1934 1990
      }
1935 1991

	
1936 1992
      /// \brief Validity checking
1937 1993
      ///
1938 1994
      /// Checks whether the iterator is invalid.
1939 1995
      bool operator==(Invalid) const { return _index == _last; }
1940 1996

	
1941 1997
      /// \brief Validity checking
1942 1998
      ///
1943 1999
      /// Checks whether the iterator is valid.
1944 2000
      bool operator!=(Invalid) const { return _index != _last; }
1945 2001

	
1946 2002
    private:
1947 2003
      const MaxWeightedMatching* _algorithm;
1948 2004
      int _last;
1949 2005
      int _index;
1950 2006
    };
1951 2007

	
1952 2008
    /// @}
1953 2009

	
1954 2010
  };
1955 2011

	
1956 2012
  /// \ingroup matching
1957 2013
  ///
1958 2014
  /// \brief Weighted perfect matching in general graphs
1959 2015
  ///
1960 2016
  /// This class provides an efficient implementation of Edmond's
1961 2017
  /// maximum weighted perfect matching algorithm. The implementation
1962 2018
  /// is based on extensive use of priority queues and provides
1963 2019
  /// \f$O(nm\log n)\f$ time complexity.
1964 2020
  ///
1965
  /// The maximum weighted matching problem is to find undirected
1966
  /// edges in the graph with maximum overall weight and no two of
1967
  /// them shares their ends and covers all nodes. The problem can be
1968
  /// formulated with the following linear program.
2021
  /// The maximum weighted perfect matching problem is to find a subset of 
2022
  /// the edges in an undirected graph with maximum overall weight for which 
2023
  /// each node has exactly one incident edge.
2024
  /// It can be formulated with the following linear program.
1969 2025
  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1970 2026
  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
1971 2027
      \quad \forall B\in\mathcal{O}\f] */
1972 2028
  /// \f[x_e \ge 0\quad \forall e\in E\f]
1973 2029
  /// \f[\max \sum_{e\in E}x_ew_e\f]
1974 2030
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1975 2031
  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
1976 2032
  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
1977 2033
  /// subsets of the nodes.
1978 2034
  ///
1979 2035
  /// The algorithm calculates an optimal matching and a proof of the
1980 2036
  /// optimality. The solution of the dual problem can be used to check
1981 2037
  /// the result of the algorithm. The dual linear problem is the
2038
  /// following.
1982 2039
  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
1983 2040
      w_{uv} \quad \forall uv\in E\f] */
1984 2041
  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1985 2042
  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1986 2043
      \frac{\vert B \vert - 1}{2}z_B\f] */
1987 2044
  ///
1988
  /// The algorithm can be executed with \c run() or the \c init() and
1989
  /// then the \c start() member functions. After it the matching can
1990
  /// be asked with \c matching() or mate() functions. The dual
1991
  /// solution can be get with \c nodeValue(), \c blossomNum() and \c
1992
  /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
1993
  /// "BlossomIt" nested class which is able to iterate on the nodes
1994
  /// of a blossom. If the value type is integral then the dual
1995
  /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
2045
  /// The algorithm can be executed with the run() function. 
2046
  /// After it the matching (the primal solution) and the dual solution
2047
  /// can be obtained using the query functions and the 
2048
  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 
2049
  /// which is able to iterate on the nodes of a blossom. 
2050
  /// If the value type is integer, then the dual solution is multiplied
2051
  /// by \ref MaxWeightedMatching::dualScale "4".
2052
  ///
2053
  /// \tparam GR The graph type the algorithm runs on.
2054
  /// \tparam WM The type edge weight map. The default type is 
2055
  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
2056
#ifdef DOXYGEN
2057
  template <typename GR, typename WM>
2058
#else
1996 2059
  template <typename GR,
1997 2060
            typename WM = typename GR::template EdgeMap<int> >
2061
#endif
1998 2062
  class MaxWeightedPerfectMatching {
1999 2063
  public:
2000 2064

	
2065
    /// The graph type of the algorithm
2001 2066
    typedef GR Graph;
2067
    /// The type of the edge weight map
2002 2068
    typedef WM WeightMap;
2069
    /// The value type of the edge weights
2003 2070
    typedef typename WeightMap::Value Value;
2004 2071

	
2005 2072
    /// \brief Scaling factor for dual solution
2006 2073
    ///
2007 2074
    /// Scaling factor for dual solution, it is equal to 4 or 1
2008 2075
    /// according to the value type.
2009 2076
    static const int dualScale =
2010 2077
      std::numeric_limits<Value>::is_integer ? 4 : 1;
2011 2078

	
2012 2079
    typedef typename Graph::template NodeMap<typename Graph::Arc>
2013 2080
    MatchingMap;
2014 2081

	
2015 2082
  private:
2016 2083

	
2017 2084
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
2018 2085

	
2019 2086
    typedef typename Graph::template NodeMap<Value> NodePotential;
2020 2087
    typedef std::vector<Node> BlossomNodeList;
2021 2088

	
2022 2089
    struct BlossomVariable {
2023 2090
      int begin, end;
2024 2091
      Value value;
2025 2092

	
2026 2093
      BlossomVariable(int _begin, int _end, Value _value)
2027 2094
        : begin(_begin), end(_end), value(_value) {}
2028 2095

	
2029 2096
    };
2030 2097

	
2031 2098
    typedef std::vector<BlossomVariable> BlossomPotential;
2032 2099

	
2033 2100
    const Graph& _graph;
2034 2101
    const WeightMap& _weight;
2035 2102

	
2036 2103
    MatchingMap* _matching;
2037 2104

	
2038 2105
    NodePotential* _node_potential;
2039 2106

	
2040 2107
    BlossomPotential _blossom_potential;
2041 2108
    BlossomNodeList _blossom_node_list;
2042 2109

	
2043 2110
    int _node_num;
2044 2111
    int _blossom_num;
2045 2112

	
2046 2113
    typedef RangeMap<int> IntIntMap;
2047 2114

	
2048 2115
    enum Status {
2049 2116
      EVEN = -1, MATCHED = 0, ODD = 1
2050 2117
    };
2051 2118

	
2052 2119
    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
2053 2120
    struct BlossomData {
2054 2121
      int tree;
2055 2122
      Status status;
2056 2123
      Arc pred, next;
2057 2124
      Value pot, offset;
2058 2125
    };
2059 2126

	
2060 2127
    IntNodeMap *_blossom_index;
2061 2128
    BlossomSet *_blossom_set;
2062 2129
    RangeMap<BlossomData>* _blossom_data;
2063 2130

	
2064 2131
    IntNodeMap *_node_index;
2065 2132
    IntArcMap *_node_heap_index;
2066 2133

	
2067 2134
    struct NodeData {
2068 2135

	
2069 2136
      NodeData(IntArcMap& node_heap_index)
2070 2137
        : heap(node_heap_index) {}
2071 2138

	
2072 2139
      int blossom;
2073 2140
      Value pot;
2074 2141
      BinHeap<Value, IntArcMap> heap;
2075 2142
      std::map<int, Arc> heap_index;
2076 2143

	
2077 2144
      int tree;
2078 2145
    };
2079 2146

	
2080 2147
    RangeMap<NodeData>* _node_data;
2081 2148

	
2082 2149
    typedef ExtendFindEnum<IntIntMap> TreeSet;
2083 2150

	
2084 2151
    IntIntMap *_tree_set_index;
2085 2152
    TreeSet *_tree_set;
2086 2153

	
2087 2154
    IntIntMap *_delta2_index;
2088 2155
    BinHeap<Value, IntIntMap> *_delta2;
2089 2156

	
2090 2157
    IntEdgeMap *_delta3_index;
2091 2158
    BinHeap<Value, IntEdgeMap> *_delta3;
2092 2159

	
2093 2160
    IntIntMap *_delta4_index;
2094 2161
    BinHeap<Value, IntIntMap> *_delta4;
2095 2162

	
2096 2163
    Value _delta_sum;
2097 2164

	
2098 2165
    void createStructures() {
2099 2166
      _node_num = countNodes(_graph);
2100 2167
      _blossom_num = _node_num * 3 / 2;
2101 2168

	
2102 2169
      if (!_matching) {
2103 2170
        _matching = new MatchingMap(_graph);
2104 2171
      }
2105 2172
      if (!_node_potential) {
2106 2173
        _node_potential = new NodePotential(_graph);
2107 2174
      }
2108 2175
      if (!_blossom_set) {
2109 2176
        _blossom_index = new IntNodeMap(_graph);
2110 2177
        _blossom_set = new BlossomSet(*_blossom_index);
2111 2178
        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2112 2179
      }
2113 2180

	
2114 2181
      if (!_node_index) {
2115 2182
        _node_index = new IntNodeMap(_graph);
2116 2183
        _node_heap_index = new IntArcMap(_graph);
2117 2184
        _node_data = new RangeMap<NodeData>(_node_num,
2118 2185
                                            NodeData(*_node_heap_index));
2119 2186
      }
2120 2187

	
2121 2188
      if (!_tree_set) {
2122 2189
        _tree_set_index = new IntIntMap(_blossom_num);
2123 2190
        _tree_set = new TreeSet(*_tree_set_index);
2124 2191
      }
2125 2192
      if (!_delta2) {
2126 2193
        _delta2_index = new IntIntMap(_blossom_num);
2127 2194
        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2128 2195
      }
2129 2196
      if (!_delta3) {
2130 2197
        _delta3_index = new IntEdgeMap(_graph);
2131 2198
        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
2132 2199
      }
2133 2200
      if (!_delta4) {
2134 2201
        _delta4_index = new IntIntMap(_blossom_num);
2135 2202
        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2136 2203
      }
2137 2204
    }
2138 2205

	
2139 2206
    void destroyStructures() {
2140 2207
      _node_num = countNodes(_graph);
2141 2208
      _blossom_num = _node_num * 3 / 2;
2142 2209

	
2143 2210
      if (_matching) {
2144 2211
        delete _matching;
2145 2212
      }
2146 2213
      if (_node_potential) {
2147 2214
        delete _node_potential;
2148 2215
      }
2149 2216
      if (_blossom_set) {
2150 2217
        delete _blossom_index;
2151 2218
        delete _blossom_set;
2152 2219
        delete _blossom_data;
2153 2220
      }
2154 2221

	
2155 2222
      if (_node_index) {
2156 2223
        delete _node_index;
2157 2224
        delete _node_heap_index;
2158 2225
        delete _node_data;
2159 2226
      }
2160 2227

	
2161 2228
      if (_tree_set) {
2162 2229
        delete _tree_set_index;
2163 2230
        delete _tree_set;
2164 2231
      }
2165 2232
      if (_delta2) {
2166 2233
        delete _delta2_index;
2167 2234
        delete _delta2;
2168 2235
      }
2169 2236
      if (_delta3) {
2170 2237
        delete _delta3_index;
2171 2238
        delete _delta3;
2172 2239
      }
2173 2240
      if (_delta4) {
2174 2241
        delete _delta4_index;
2175 2242
        delete _delta4;
2176 2243
      }
2177 2244
    }
2178 2245

	
2179 2246
    void matchedToEven(int blossom, int tree) {
2180 2247
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2181 2248
        _delta2->erase(blossom);
2182 2249
      }
2183 2250

	
2184 2251
      if (!_blossom_set->trivial(blossom)) {
2185 2252
        (*_blossom_data)[blossom].pot -=
2186 2253
          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2187 2254
      }
2188 2255

	
2189 2256
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2190 2257
           n != INVALID; ++n) {
2191 2258

	
2192 2259
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
2193 2260
        int ni = (*_node_index)[n];
2194 2261

	
... ...
@@ -2629,479 +2696,512 @@
2629 2696
      Arc next = (*_blossom_data)[blossom].next;
2630 2697
      Arc pred = (*_blossom_data)[blossom].pred;
2631 2698

	
2632 2699
      int tree = _tree_set->find(blossom);
2633 2700

	
2634 2701
      (*_blossom_data)[blossom].status = MATCHED;
2635 2702
      oddToMatched(blossom);
2636 2703
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2637 2704
        _delta2->erase(blossom);
2638 2705
      }
2639 2706

	
2640 2707
      std::vector<int> subblossoms;
2641 2708
      _blossom_set->split(blossom, std::back_inserter(subblossoms));
2642 2709

	
2643 2710
      Value offset = (*_blossom_data)[blossom].offset;
2644 2711
      int b = _blossom_set->find(_graph.source(pred));
2645 2712
      int d = _blossom_set->find(_graph.source(next));
2646 2713

	
2647 2714
      int ib = -1, id = -1;
2648 2715
      for (int i = 0; i < int(subblossoms.size()); ++i) {
2649 2716
        if (subblossoms[i] == b) ib = i;
2650 2717
        if (subblossoms[i] == d) id = i;
2651 2718

	
2652 2719
        (*_blossom_data)[subblossoms[i]].offset = offset;
2653 2720
        if (!_blossom_set->trivial(subblossoms[i])) {
2654 2721
          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2655 2722
        }
2656 2723
        if (_blossom_set->classPrio(subblossoms[i]) !=
2657 2724
            std::numeric_limits<Value>::max()) {
2658 2725
          _delta2->push(subblossoms[i],
2659 2726
                        _blossom_set->classPrio(subblossoms[i]) -
2660 2727
                        (*_blossom_data)[subblossoms[i]].offset);
2661 2728
        }
2662 2729
      }
2663 2730

	
2664 2731
      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2665 2732
        for (int i = (id + 1) % subblossoms.size();
2666 2733
             i != ib; i = (i + 2) % subblossoms.size()) {
2667 2734
          int sb = subblossoms[i];
2668 2735
          int tb = subblossoms[(i + 1) % subblossoms.size()];
2669 2736
          (*_blossom_data)[sb].next =
2670 2737
            _graph.oppositeArc((*_blossom_data)[tb].next);
2671 2738
        }
2672 2739

	
2673 2740
        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2674 2741
          int sb = subblossoms[i];
2675 2742
          int tb = subblossoms[(i + 1) % subblossoms.size()];
2676 2743
          int ub = subblossoms[(i + 2) % subblossoms.size()];
2677 2744

	
2678 2745
          (*_blossom_data)[sb].status = ODD;
2679 2746
          matchedToOdd(sb);
2680 2747
          _tree_set->insert(sb, tree);
2681 2748
          (*_blossom_data)[sb].pred = pred;
2682 2749
          (*_blossom_data)[sb].next =
2683 2750
                           _graph.oppositeArc((*_blossom_data)[tb].next);
2684 2751

	
2685 2752
          pred = (*_blossom_data)[ub].next;
2686 2753

	
2687 2754
          (*_blossom_data)[tb].status = EVEN;
2688 2755
          matchedToEven(tb, tree);
2689 2756
          _tree_set->insert(tb, tree);
2690 2757
          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2691 2758
        }
2692 2759

	
2693 2760
        (*_blossom_data)[subblossoms[id]].status = ODD;
2694 2761
        matchedToOdd(subblossoms[id]);
2695 2762
        _tree_set->insert(subblossoms[id], tree);
2696 2763
        (*_blossom_data)[subblossoms[id]].next = next;
2697 2764
        (*_blossom_data)[subblossoms[id]].pred = pred;
2698 2765

	
2699 2766
      } else {
2700 2767

	
2701 2768
        for (int i = (ib + 1) % subblossoms.size();
2702 2769
             i != id; i = (i + 2) % subblossoms.size()) {
2703 2770
          int sb = subblossoms[i];
2704 2771
          int tb = subblossoms[(i + 1) % subblossoms.size()];
2705 2772
          (*_blossom_data)[sb].next =
2706 2773
            _graph.oppositeArc((*_blossom_data)[tb].next);
2707 2774
        }
2708 2775

	
2709 2776
        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2710 2777
          int sb = subblossoms[i];
2711 2778
          int tb = subblossoms[(i + 1) % subblossoms.size()];
2712 2779
          int ub = subblossoms[(i + 2) % subblossoms.size()];
2713 2780

	
2714 2781
          (*_blossom_data)[sb].status = ODD;
2715 2782
          matchedToOdd(sb);
2716 2783
          _tree_set->insert(sb, tree);
2717 2784
          (*_blossom_data)[sb].next = next;
2718 2785
          (*_blossom_data)[sb].pred =
2719 2786
            _graph.oppositeArc((*_blossom_data)[tb].next);
2720 2787

	
2721 2788
          (*_blossom_data)[tb].status = EVEN;
2722 2789
          matchedToEven(tb, tree);
2723 2790
          _tree_set->insert(tb, tree);
2724 2791
          (*_blossom_data)[tb].pred =
2725 2792
            (*_blossom_data)[tb].next =
2726 2793
            _graph.oppositeArc((*_blossom_data)[ub].next);
2727 2794
          next = (*_blossom_data)[ub].next;
2728 2795
        }
2729 2796

	
2730 2797
        (*_blossom_data)[subblossoms[ib]].status = ODD;
2731 2798
        matchedToOdd(subblossoms[ib]);
2732 2799
        _tree_set->insert(subblossoms[ib], tree);
2733 2800
        (*_blossom_data)[subblossoms[ib]].next = next;
2734 2801
        (*_blossom_data)[subblossoms[ib]].pred = pred;
2735 2802
      }
2736 2803
      _tree_set->erase(blossom);
2737 2804
    }
2738 2805

	
2739 2806
    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2740 2807
      if (_blossom_set->trivial(blossom)) {
2741 2808
        int bi = (*_node_index)[base];
2742 2809
        Value pot = (*_node_data)[bi].pot;
2743 2810

	
2744 2811
        (*_matching)[base] = matching;
2745 2812
        _blossom_node_list.push_back(base);
2746 2813
        (*_node_potential)[base] = pot;
2747 2814
      } else {
2748 2815

	
2749 2816
        Value pot = (*_blossom_data)[blossom].pot;
2750 2817
        int bn = _blossom_node_list.size();
2751 2818

	
2752 2819
        std::vector<int> subblossoms;
2753 2820
        _blossom_set->split(blossom, std::back_inserter(subblossoms));
2754 2821
        int b = _blossom_set->find(base);
2755 2822
        int ib = -1;
2756 2823
        for (int i = 0; i < int(subblossoms.size()); ++i) {
2757 2824
          if (subblossoms[i] == b) { ib = i; break; }
2758 2825
        }
2759 2826

	
2760 2827
        for (int i = 1; i < int(subblossoms.size()); i += 2) {
2761 2828
          int sb = subblossoms[(ib + i) % subblossoms.size()];
2762 2829
          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2763 2830

	
2764 2831
          Arc m = (*_blossom_data)[tb].next;
2765 2832
          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2766 2833
          extractBlossom(tb, _graph.source(m), m);
2767 2834
        }
2768 2835
        extractBlossom(subblossoms[ib], base, matching);
2769 2836

	
2770 2837
        int en = _blossom_node_list.size();
2771 2838

	
2772 2839
        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2773 2840
      }
2774 2841
    }
2775 2842

	
2776 2843
    void extractMatching() {
2777 2844
      std::vector<int> blossoms;
2778 2845
      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2779 2846
        blossoms.push_back(c);
2780 2847
      }
2781 2848

	
2782 2849
      for (int i = 0; i < int(blossoms.size()); ++i) {
2783 2850

	
2784 2851
        Value offset = (*_blossom_data)[blossoms[i]].offset;
2785 2852
        (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2786 2853
        for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2787 2854
             n != INVALID; ++n) {
2788 2855
          (*_node_data)[(*_node_index)[n]].pot -= offset;
2789 2856
        }
2790 2857

	
2791 2858
        Arc matching = (*_blossom_data)[blossoms[i]].next;
2792 2859
        Node base = _graph.source(matching);
2793 2860
        extractBlossom(blossoms[i], base, matching);
2794 2861
      }
2795 2862
    }
2796 2863

	
2797 2864
  public:
2798 2865

	
2799 2866
    /// \brief Constructor
2800 2867
    ///
2801 2868
    /// Constructor.
2802 2869
    MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2803 2870
      : _graph(graph), _weight(weight), _matching(0),
2804 2871
        _node_potential(0), _blossom_potential(), _blossom_node_list(),
2805 2872
        _node_num(0), _blossom_num(0),
2806 2873

	
2807 2874
        _blossom_index(0), _blossom_set(0), _blossom_data(0),
2808 2875
        _node_index(0), _node_heap_index(0), _node_data(0),
2809 2876
        _tree_set_index(0), _tree_set(0),
2810 2877

	
2811 2878
        _delta2_index(0), _delta2(0),
2812 2879
        _delta3_index(0), _delta3(0),
2813 2880
        _delta4_index(0), _delta4(0),
2814 2881

	
2815 2882
        _delta_sum() {}
2816 2883

	
2817 2884
    ~MaxWeightedPerfectMatching() {
2818 2885
      destroyStructures();
2819 2886
    }
2820 2887

	
2821
    /// \name Execution control
2888
    /// \name Execution Control
2822 2889
    /// The simplest way to execute the algorithm is to use the
2823
    /// \c run() member function.
2890
    /// \ref run() member function.
2824 2891

	
2825 2892
    ///@{
2826 2893

	
2827 2894
    /// \brief Initialize the algorithm
2828 2895
    ///
2829
    /// Initialize the algorithm
2896
    /// This function initializes the algorithm.
2830 2897
    void init() {
2831 2898
      createStructures();
2832 2899

	
2833 2900
      for (ArcIt e(_graph); e != INVALID; ++e) {
2834 2901
        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
2835 2902
      }
2836 2903
      for (EdgeIt e(_graph); e != INVALID; ++e) {
2837 2904
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
2838 2905
      }
2839 2906
      for (int i = 0; i < _blossom_num; ++i) {
2840 2907
        (*_delta2_index)[i] = _delta2->PRE_HEAP;
2841 2908
        (*_delta4_index)[i] = _delta4->PRE_HEAP;
2842 2909
      }
2843 2910

	
2844 2911
      int index = 0;
2845 2912
      for (NodeIt n(_graph); n != INVALID; ++n) {
2846 2913
        Value max = - std::numeric_limits<Value>::max();
2847 2914
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2848 2915
          if (_graph.target(e) == n) continue;
2849 2916
          if ((dualScale * _weight[e]) / 2 > max) {
2850 2917
            max = (dualScale * _weight[e]) / 2;
2851 2918
          }
2852 2919
        }
2853 2920
        (*_node_index)[n] = index;
2854 2921
        (*_node_data)[index].pot = max;
2855 2922
        int blossom =
2856 2923
          _blossom_set->insert(n, std::numeric_limits<Value>::max());
2857 2924

	
2858 2925
        _tree_set->insert(blossom);
2859 2926

	
2860 2927
        (*_blossom_data)[blossom].status = EVEN;
2861 2928
        (*_blossom_data)[blossom].pred = INVALID;
2862 2929
        (*_blossom_data)[blossom].next = INVALID;
2863 2930
        (*_blossom_data)[blossom].pot = 0;
2864 2931
        (*_blossom_data)[blossom].offset = 0;
2865 2932
        ++index;
2866 2933
      }
2867 2934
      for (EdgeIt e(_graph); e != INVALID; ++e) {
2868 2935
        int si = (*_node_index)[_graph.u(e)];
2869 2936
        int ti = (*_node_index)[_graph.v(e)];
2870 2937
        if (_graph.u(e) != _graph.v(e)) {
2871 2938
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
2872 2939
                            dualScale * _weight[e]) / 2);
2873 2940
        }
2874 2941
      }
2875 2942
    }
2876 2943

	
2877
    /// \brief Starts the algorithm
2944
    /// \brief Start the algorithm
2878 2945
    ///
2879
    /// Starts the algorithm
2946
    /// This function starts the algorithm.
2947
    ///
2948
    /// \pre \ref init() must be called before using this function.
2880 2949
    bool start() {
2881 2950
      enum OpType {
2882 2951
        D2, D3, D4
2883 2952
      };
2884 2953

	
2885 2954
      int unmatched = _node_num;
2886 2955
      while (unmatched > 0) {
2887 2956
        Value d2 = !_delta2->empty() ?
2888 2957
          _delta2->prio() : std::numeric_limits<Value>::max();
2889 2958

	
2890 2959
        Value d3 = !_delta3->empty() ?
2891 2960
          _delta3->prio() : std::numeric_limits<Value>::max();
2892 2961

	
2893 2962
        Value d4 = !_delta4->empty() ?
2894 2963
          _delta4->prio() : std::numeric_limits<Value>::max();
2895 2964

	
2896 2965
        _delta_sum = d2; OpType ot = D2;
2897 2966
        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
2898 2967
        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
2899 2968

	
2900 2969
        if (_delta_sum == std::numeric_limits<Value>::max()) {
2901 2970
          return false;
2902 2971
        }
2903 2972

	
2904 2973
        switch (ot) {
2905 2974
        case D2:
2906 2975
          {
2907 2976
            int blossom = _delta2->top();
2908 2977
            Node n = _blossom_set->classTop(blossom);
2909 2978
            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
2910 2979
            extendOnArc(e);
2911 2980
          }
2912 2981
          break;
2913 2982
        case D3:
2914 2983
          {
2915 2984
            Edge e = _delta3->top();
2916 2985

	
2917 2986
            int left_blossom = _blossom_set->find(_graph.u(e));
2918 2987
            int right_blossom = _blossom_set->find(_graph.v(e));
2919 2988

	
2920 2989
            if (left_blossom == right_blossom) {
2921 2990
              _delta3->pop();
2922 2991
            } else {
2923 2992
              int left_tree = _tree_set->find(left_blossom);
2924 2993
              int right_tree = _tree_set->find(right_blossom);
2925 2994

	
2926 2995
              if (left_tree == right_tree) {
2927 2996
                shrinkOnEdge(e, left_tree);
2928 2997
              } else {
2929 2998
                augmentOnEdge(e);
2930 2999
                unmatched -= 2;
2931 3000
              }
2932 3001
            }
2933 3002
          } break;
2934 3003
        case D4:
2935 3004
          splitBlossom(_delta4->top());
2936 3005
          break;
2937 3006
        }
2938 3007
      }
2939 3008
      extractMatching();
2940 3009
      return true;
2941 3010
    }
2942 3011

	
2943
    /// \brief Runs %MaxWeightedPerfectMatching algorithm.
3012
    /// \brief Run the algorithm.
2944 3013
    ///
2945
    /// This method runs the %MaxWeightedPerfectMatching algorithm.
3014
    /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
2946 3015
    ///
2947
    /// \note mwm.run() is just a shortcut of the following code.
3016
    /// \note mwpm.run() is just a shortcut of the following code.
2948 3017
    /// \code
2949
    ///   mwm.init();
2950
    ///   mwm.start();
3018
    ///   mwpm.init();
3019
    ///   mwpm.start();
2951 3020
    /// \endcode
2952 3021
    bool run() {
2953 3022
      init();
2954 3023
      return start();
2955 3024
    }
2956 3025

	
2957 3026
    /// @}
2958 3027

	
2959
    /// \name Primal solution
2960
    /// Functions to get the primal solution, ie. the matching.
3028
    /// \name Primal Solution
3029
    /// Functions to get the primal solution, i.e. the maximum weighted 
3030
    /// perfect matching.\n
3031
    /// Either \ref run() or \ref start() function should be called before
3032
    /// using them.
2961 3033

	
2962 3034
    /// @{
2963 3035

	
2964
    /// \brief Returns the matching value.
3036
    /// \brief Return the weight of the matching.
2965 3037
    ///
2966
    /// Returns the matching value.
3038
    /// This function returns the weight of the found matching.
3039
    ///
3040
    /// \pre Either run() or start() must be called before using this function.
2967 3041
    Value matchingValue() const {
2968 3042
      Value sum = 0;
2969 3043
      for (NodeIt n(_graph); n != INVALID; ++n) {
2970 3044
        if ((*_matching)[n] != INVALID) {
2971 3045
          sum += _weight[(*_matching)[n]];
2972 3046
        }
2973 3047
      }
2974 3048
      return sum /= 2;
2975 3049
    }
2976 3050

	
2977
    /// \brief Returns true when the edge is in the matching.
3051
    /// \brief Return \c true if the given edge is in the matching.
2978 3052
    ///
2979
    /// Returns true when the edge is in the matching.
3053
    /// This function returns \c true if the given edge is in the found 
3054
    /// matching.
3055
    ///
3056
    /// \pre Either run() or start() must be called before using this function.
2980 3057
    bool matching(const Edge& edge) const {
2981 3058
      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
2982 3059
    }
2983 3060

	
2984
    /// \brief Returns the incident matching edge.
3061
    /// \brief Return the matching arc (or edge) incident to the given node.
2985 3062
    ///
2986
    /// Returns the incident matching arc from given edge.
3063
    /// This function returns the matching arc (or edge) incident to the
3064
    /// given node in the found matching or \c INVALID if the node is 
3065
    /// not covered by the matching.
3066
    ///
3067
    /// \pre Either run() or start() must be called before using this function.
2987 3068
    Arc matching(const Node& node) const {
2988 3069
      return (*_matching)[node];
2989 3070
    }
2990 3071

	
2991
    /// \brief Returns the mate of the node.
3072
    /// \brief Return the mate of the given node.
2992 3073
    ///
2993
    /// Returns the adjancent node in a mathcing arc.
3074
    /// This function returns the mate of the given node in the found 
3075
    /// matching or \c INVALID if the node is not covered by the matching.
3076
    ///
3077
    /// \pre Either run() or start() must be called before using this function.
2994 3078
    Node mate(const Node& node) const {
2995 3079
      return _graph.target((*_matching)[node]);
2996 3080
    }
2997 3081

	
2998 3082
    /// @}
2999 3083

	
3000
    /// \name Dual solution
3001
    /// Functions to get the dual solution.
3084
    /// \name Dual Solution
3085
    /// Functions to get the dual solution.\n
3086
    /// Either \ref run() or \ref start() function should be called before
3087
    /// using them.
3002 3088

	
3003 3089
    /// @{
3004 3090

	
3005
    /// \brief Returns the value of the dual solution.
3091
    /// \brief Return the value of the dual solution.
3006 3092
    ///
3007
    /// Returns the value of the dual solution. It should be equal to
3008
    /// the primal value scaled by \ref dualScale "dual scale".
3093
    /// This function returns the value of the dual solution. 
3094
    /// It should be equal to the primal value scaled by \ref dualScale 
3095
    /// "dual scale".
3096
    ///
3097
    /// \pre Either run() or start() must be called before using this function.
3009 3098
    Value dualValue() const {
3010 3099
      Value sum = 0;
3011 3100
      for (NodeIt n(_graph); n != INVALID; ++n) {
3012 3101
        sum += nodeValue(n);
3013 3102
      }
3014 3103
      for (int i = 0; i < blossomNum(); ++i) {
3015 3104
        sum += blossomValue(i) * (blossomSize(i) / 2);
3016 3105
      }
3017 3106
      return sum;
3018 3107
    }
3019 3108

	
3020
    /// \brief Returns the value of the node.
3109
    /// \brief Return the dual value (potential) of the given node.
3021 3110
    ///
3022
    /// Returns the the value of the node.
3111
    /// This function returns the dual value (potential) of the given node.
3112
    ///
3113
    /// \pre Either run() or start() must be called before using this function.
3023 3114
    Value nodeValue(const Node& n) const {
3024 3115
      return (*_node_potential)[n];
3025 3116
    }
3026 3117

	
3027
    /// \brief Returns the number of the blossoms in the basis.
3118
    /// \brief Return the number of the blossoms in the basis.
3028 3119
    ///
3029
    /// Returns the number of the blossoms in the basis.
3120
    /// This function returns the number of the blossoms in the basis.
3121
    ///
3122
    /// \pre Either run() or start() must be called before using this function.
3030 3123
    /// \see BlossomIt
3031 3124
    int blossomNum() const {
3032 3125
      return _blossom_potential.size();
3033 3126
    }
3034 3127

	
3035

	
3036
    /// \brief Returns the number of the nodes in the blossom.
3128
    /// \brief Return the number of the nodes in the given blossom.
3037 3129
    ///
3038
    /// Returns the number of the nodes in the blossom.
3130
    /// This function returns the number of the nodes in the given blossom.
3131
    ///
3132
    /// \pre Either run() or start() must be called before using this function.
3133
    /// \see BlossomIt
3039 3134
    int blossomSize(int k) const {
3040 3135
      return _blossom_potential[k].end - _blossom_potential[k].begin;
3041 3136
    }
3042 3137

	
3043
    /// \brief Returns the value of the blossom.
3138
    /// \brief Return the dual value (ptential) of the given blossom.
3044 3139
    ///
3045
    /// Returns the the value of the blossom.
3046
    /// \see BlossomIt
3140
    /// This function returns the dual value (ptential) of the given blossom.
3141
    ///
3142
    /// \pre Either run() or start() must be called before using this function.
3047 3143
    Value blossomValue(int k) const {
3048 3144
      return _blossom_potential[k].value;
3049 3145
    }
3050 3146

	
3051
    /// \brief Iterator for obtaining the nodes of the blossom.
3147
    /// \brief Iterator for obtaining the nodes of a blossom.
3052 3148
    ///
3053
    /// Iterator for obtaining the nodes of the blossom. This class
3054
    /// provides a common lemon style iterator for listing a
3055
    /// subset of the nodes.
3149
    /// This class provides an iterator for obtaining the nodes of the 
3150
    /// given blossom. It lists a subset of the nodes.
3151
    /// Before using this iterator, you must allocate a 
3152
    /// MaxWeightedPerfectMatching class and execute it.
3056 3153
    class BlossomIt {
3057 3154
    public:
3058 3155

	
3059 3156
      /// \brief Constructor.
3060 3157
      ///
3061
      /// Constructor to get the nodes of the variable.
3158
      /// Constructor to get the nodes of the given variable.
3159
      ///
3160
      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 
3161
      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 
3162
      /// must be called before initializing this iterator.
3062 3163
      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
3063 3164
        : _algorithm(&algorithm)
3064 3165
      {
3065 3166
        _index = _algorithm->_blossom_potential[variable].begin;
3066 3167
        _last = _algorithm->_blossom_potential[variable].end;
3067 3168
      }
3068 3169

	
3069
      /// \brief Conversion to node.
3170
      /// \brief Conversion to \c Node.
3070 3171
      ///
3071
      /// Conversion to node.
3172
      /// Conversion to \c Node.
3072 3173
      operator Node() const {
3073 3174
        return _algorithm->_blossom_node_list[_index];
3074 3175
      }
3075 3176

	
3076 3177
      /// \brief Increment operator.
3077 3178
      ///
3078 3179
      /// Increment operator.
3079 3180
      BlossomIt& operator++() {
3080 3181
        ++_index;
3081 3182
        return *this;
3082 3183
      }
3083 3184

	
3084 3185
      /// \brief Validity checking
3085 3186
      ///
3086
      /// Checks whether the iterator is invalid.
3187
      /// This function checks whether the iterator is invalid.
3087 3188
      bool operator==(Invalid) const { return _index == _last; }
3088 3189

	
3089 3190
      /// \brief Validity checking
3090 3191
      ///
3091
      /// Checks whether the iterator is valid.
3192
      /// This function checks whether the iterator is valid.
3092 3193
      bool operator!=(Invalid) const { return _index != _last; }
3093 3194

	
3094 3195
    private:
3095 3196
      const MaxWeightedPerfectMatching* _algorithm;
3096 3197
      int _last;
3097 3198
      int _index;
3098 3199
    };
3099 3200

	
3100 3201
    /// @}
3101 3202

	
3102 3203
  };
3103 3204

	
3104

	
3105 3205
} //END OF NAMESPACE LEMON
3106 3206

	
3107 3207
#endif //LEMON_MAX_MATCHING_H
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#include <iostream>
20 20
#include <sstream>
21 21
#include <vector>
22 22
#include <queue>
23
#include <lemon/math.h>
24 23
#include <cstdlib>
25 24

	
26 25
#include <lemon/max_matching.h>
27 26
#include <lemon/smart_graph.h>
27
#include <lemon/concepts/graph.h>
28
#include <lemon/concepts/maps.h>
28 29
#include <lemon/lgf_reader.h>
30
#include <lemon/math.h>
29 31

	
30 32
#include "test_tools.h"
31 33

	
32 34
using namespace std;
33 35
using namespace lemon;
34 36

	
35 37
GRAPH_TYPEDEFS(SmartGraph);
36 38

	
37 39

	
38 40
const int lgfn = 3;
39 41
const std::string lgf[lgfn] = {
40 42
  "@nodes\n"
41 43
  "label\n"
42 44
  "0\n"
43 45
  "1\n"
44 46
  "2\n"
45 47
  "3\n"
46 48
  "4\n"
47 49
  "5\n"
48 50
  "6\n"
49 51
  "7\n"
50 52
  "@edges\n"
51 53
  "     label  weight\n"
52 54
  "7 4  0      984\n"
53 55
  "0 7  1      73\n"
54 56
  "7 1  2      204\n"
55 57
  "2 3  3      583\n"
56 58
  "2 7  4      565\n"
57 59
  "2 1  5      582\n"
58 60
  "0 4  6      551\n"
59 61
  "2 5  7      385\n"
60 62
  "1 5  8      561\n"
61 63
  "5 3  9      484\n"
62 64
  "7 5  10     904\n"
63 65
  "3 6  11     47\n"
64 66
  "7 6  12     888\n"
65 67
  "3 0  13     747\n"
66 68
  "6 1  14     310\n",
67 69

	
68 70
  "@nodes\n"
69 71
  "label\n"
70 72
  "0\n"
71 73
  "1\n"
72 74
  "2\n"
73 75
  "3\n"
74 76
  "4\n"
75 77
  "5\n"
76 78
  "6\n"
77 79
  "7\n"
78 80
  "@edges\n"
79 81
  "     label  weight\n"
80 82
  "2 5  0      710\n"
81 83
  "0 5  1      241\n"
82 84
  "2 4  2      856\n"
83 85
  "2 6  3      762\n"
84 86
  "4 1  4      747\n"
85 87
  "6 1  5      962\n"
86 88
  "4 7  6      723\n"
87 89
  "1 7  7      661\n"
88 90
  "2 3  8      376\n"
89 91
  "1 0  9      416\n"
90 92
  "6 7  10     391\n",
91 93

	
92 94
  "@nodes\n"
93 95
  "label\n"
94 96
  "0\n"
95 97
  "1\n"
96 98
  "2\n"
97 99
  "3\n"
98 100
  "4\n"
99 101
  "5\n"
100 102
  "6\n"
101 103
  "7\n"
102 104
  "@edges\n"
103 105
  "     label  weight\n"
104 106
  "6 2  0      553\n"
105 107
  "0 7  1      653\n"
106 108
  "6 3  2      22\n"
107 109
  "4 7  3      846\n"
108 110
  "7 2  4      981\n"
109 111
  "7 6  5      250\n"
110 112
  "5 2  6      539\n",
111 113
};
112 114

	
115
void checkMaxMatchingCompile()
116
{
117
  typedef concepts::Graph Graph;
118
  typedef Graph::Node Node;
119
  typedef Graph::Edge Edge;
120
  typedef Graph::EdgeMap<bool> MatMap;
121

	
122
  Graph g;
123
  Node n;
124
  Edge e;
125
  MatMap mat(g);
126

	
127
  MaxMatching<Graph> mat_test(g);
128
  const MaxMatching<Graph>&
129
    const_mat_test = mat_test;
130

	
131
  mat_test.init();
132
  mat_test.greedyInit();
133
  mat_test.matchingInit(mat);
134
  mat_test.startSparse();
135
  mat_test.startDense();
136
  mat_test.run();
137
  
138
  const_mat_test.matchingSize();
139
  const_mat_test.matching(e);
140
  const_mat_test.matching(n);
141
  const_mat_test.mate(n);
142

	
143
  MaxMatching<Graph>::Status stat = 
144
    const_mat_test.decomposition(n);
145
  const_mat_test.barrier(n);
146
  
147
  ignore_unused_variable_warning(stat);
148
}
149

	
150
void checkMaxWeightedMatchingCompile()
151
{
152
  typedef concepts::Graph Graph;
153
  typedef Graph::Node Node;
154
  typedef Graph::Edge Edge;
155
  typedef Graph::EdgeMap<int> WeightMap;
156

	
157
  Graph g;
158
  Node n;
159
  Edge e;
160
  WeightMap w(g);
161

	
162
  MaxWeightedMatching<Graph> mat_test(g, w);
163
  const MaxWeightedMatching<Graph>&
164
    const_mat_test = mat_test;
165

	
166
  mat_test.init();
167
  mat_test.start();
168
  mat_test.run();
169
  
170
  const_mat_test.matchingValue();
171
  const_mat_test.matchingSize();
172
  const_mat_test.matching(e);
173
  const_mat_test.matching(n);
174
  const_mat_test.mate(n);
175
  
176
  int k = 0;
177
  const_mat_test.dualValue();
178
  const_mat_test.nodeValue(n);
179
  const_mat_test.blossomNum();
180
  const_mat_test.blossomSize(k);
181
  const_mat_test.blossomValue(k);
182
}
183

	
184
void checkMaxWeightedPerfectMatchingCompile()
185
{
186
  typedef concepts::Graph Graph;
187
  typedef Graph::Node Node;
188
  typedef Graph::Edge Edge;
189
  typedef Graph::EdgeMap<int> WeightMap;
190

	
191
  Graph g;
192
  Node n;
193
  Edge e;
194
  WeightMap w(g);
195

	
196
  MaxWeightedPerfectMatching<Graph> mat_test(g, w);
197
  const MaxWeightedPerfectMatching<Graph>&
198
    const_mat_test = mat_test;
199

	
200
  mat_test.init();
201
  mat_test.start();
202
  mat_test.run();
203
  
204
  const_mat_test.matchingValue();
205
  const_mat_test.matching(e);
206
  const_mat_test.matching(n);
207
  const_mat_test.mate(n);
208
  
209
  int k = 0;
210
  const_mat_test.dualValue();
211
  const_mat_test.nodeValue(n);
212
  const_mat_test.blossomNum();
213
  const_mat_test.blossomSize(k);
214
  const_mat_test.blossomValue(k);
215
}
216

	
113 217
void checkMatching(const SmartGraph& graph,
114 218
                   const MaxMatching<SmartGraph>& mm) {
115 219
  int num = 0;
116 220

	
117 221
  IntNodeMap comp_index(graph);
118 222
  UnionFind<IntNodeMap> comp(comp_index);
119 223

	
120 224
  int barrier_num = 0;
121 225

	
122 226
  for (NodeIt n(graph); n != INVALID; ++n) {
123 227
    check(mm.decomposition(n) == MaxMatching<SmartGraph>::EVEN ||
124 228
          mm.matching(n) != INVALID, "Wrong Gallai-Edmonds decomposition");
125 229
    if (mm.decomposition(n) == MaxMatching<SmartGraph>::ODD) {
126 230
      ++barrier_num;
127 231
    } else {
128 232
      comp.insert(n);
129 233
    }
130 234
  }
131 235

	
132 236
  for (EdgeIt e(graph); e != INVALID; ++e) {
133 237
    if (mm.matching(e)) {
134 238
      check(e == mm.matching(graph.u(e)), "Wrong matching");
135 239
      check(e == mm.matching(graph.v(e)), "Wrong matching");
136 240
      ++num;
137 241
    }
138 242
    check(mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::EVEN ||
139 243
          mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::MATCHED,
140 244
          "Wrong Gallai-Edmonds decomposition");
141 245

	
142 246
    check(mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::EVEN ||
143 247
          mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::MATCHED,
144 248
          "Wrong Gallai-Edmonds decomposition");
145 249

	
146 250
    if (mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::ODD &&
147 251
        mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::ODD) {
148 252
      comp.join(graph.u(e), graph.v(e));
149 253
    }
150 254
  }
151 255

	
152 256
  std::set<int> comp_root;
153 257
  int odd_comp_num = 0;
154 258
  for (NodeIt n(graph); n != INVALID; ++n) {
155 259
    if (mm.decomposition(n) != MaxMatching<SmartGraph>::ODD) {
156 260
      int root = comp.find(n);
157 261
      if (comp_root.find(root) == comp_root.end()) {
158 262
        comp_root.insert(root);
159 263
        if (comp.size(n) % 2 == 1) {
160 264
          ++odd_comp_num;
161 265
        }
162 266
      }
163 267
    }
164 268
  }
165 269

	
166 270
  check(mm.matchingSize() == num, "Wrong matching");
167 271
  check(2 * num == countNodes(graph) - (odd_comp_num - barrier_num),
168 272
         "Wrong matching");
169 273
  return;
170 274
}
171 275

	
172 276
void checkWeightedMatching(const SmartGraph& graph,
173 277
                   const SmartGraph::EdgeMap<int>& weight,
174 278
                   const MaxWeightedMatching<SmartGraph>& mwm) {
175 279
  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
176 280
    if (graph.u(e) == graph.v(e)) continue;
177 281
    int rw = mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
178 282

	
179 283
    for (int i = 0; i < mwm.blossomNum(); ++i) {
180 284
      bool s = false, t = false;
181 285
      for (MaxWeightedMatching<SmartGraph>::BlossomIt n(mwm, i);
182 286
           n != INVALID; ++n) {
183 287
        if (graph.u(e) == n) s = true;
184 288
        if (graph.v(e) == n) t = true;
185 289
      }
186 290
      if (s == true && t == true) {
187 291
        rw += mwm.blossomValue(i);
188 292
      }
189 293
    }
190 294
    rw -= weight[e] * mwm.dualScale;
191 295

	
192 296
    check(rw >= 0, "Negative reduced weight");
193 297
    check(rw == 0 || !mwm.matching(e),
194 298
          "Non-zero reduced weight on matching edge");
195 299
  }
196 300

	
197 301
  int pv = 0;
198 302
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
199 303
    if (mwm.matching(n) != INVALID) {
200 304
      check(mwm.nodeValue(n) >= 0, "Invalid node value");
201 305
      pv += weight[mwm.matching(n)];
202 306
      SmartGraph::Node o = graph.target(mwm.matching(n));
203 307
      check(mwm.mate(n) == o, "Invalid matching");
204 308
      check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
205 309
            "Invalid matching");
206 310
    } else {
207 311
      check(mwm.mate(n) == INVALID, "Invalid matching");
208 312
      check(mwm.nodeValue(n) == 0, "Invalid matching");
209 313
    }
210 314
  }
211 315

	
212 316
  int dv = 0;
213 317
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
214 318
    dv += mwm.nodeValue(n);
215 319
  }
216 320

	
217 321
  for (int i = 0; i < mwm.blossomNum(); ++i) {
218 322
    check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
219 323
    check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
220 324
    dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
221 325
  }
222 326

	
223 327
  check(pv * mwm.dualScale == dv * 2, "Wrong duality");
224 328

	
225 329
  return;
226 330
}
227 331

	
228 332
void checkWeightedPerfectMatching(const SmartGraph& graph,
229 333
                          const SmartGraph::EdgeMap<int>& weight,
230 334
                          const MaxWeightedPerfectMatching<SmartGraph>& mwpm) {
231 335
  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
232 336
    if (graph.u(e) == graph.v(e)) continue;
233 337
    int rw = mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
234 338

	
235 339
    for (int i = 0; i < mwpm.blossomNum(); ++i) {
236 340
      bool s = false, t = false;
237 341
      for (MaxWeightedPerfectMatching<SmartGraph>::BlossomIt n(mwpm, i);
238 342
           n != INVALID; ++n) {
239 343
        if (graph.u(e) == n) s = true;
240 344
        if (graph.v(e) == n) t = true;
241 345
      }
242 346
      if (s == true && t == true) {
243 347
        rw += mwpm.blossomValue(i);
244 348
      }
245 349
    }
246 350
    rw -= weight[e] * mwpm.dualScale;
247 351

	
248 352
    check(rw >= 0, "Negative reduced weight");
249 353
    check(rw == 0 || !mwpm.matching(e),
250 354
          "Non-zero reduced weight on matching edge");
251 355
  }
252 356

	
253 357
  int pv = 0;
254 358
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
255 359
    check(mwpm.matching(n) != INVALID, "Non perfect");
256 360
    pv += weight[mwpm.matching(n)];
257 361
    SmartGraph::Node o = graph.target(mwpm.matching(n));
258 362
    check(mwpm.mate(n) == o, "Invalid matching");
259 363
    check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
260 364
          "Invalid matching");
261 365
  }
262 366

	
263 367
  int dv = 0;
264 368
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
265 369
    dv += mwpm.nodeValue(n);
266 370
  }
267 371

	
268 372
  for (int i = 0; i < mwpm.blossomNum(); ++i) {
269 373
    check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
270 374
    check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
271 375
    dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
272 376
  }
273 377

	
274 378
  check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
275 379

	
276 380
  return;
277 381
}
278 382

	
279 383

	
280 384
int main() {
281 385

	
282 386
  for (int i = 0; i < lgfn; ++i) {
283 387
    SmartGraph graph;
284 388
    SmartGraph::EdgeMap<int> weight(graph);
285 389

	
286 390
    istringstream lgfs(lgf[i]);
287 391
    graphReader(graph, lgfs).
288 392
      edgeMap("weight", weight).run();
289 393

	
290 394
    MaxMatching<SmartGraph> mm(graph);
291 395
    mm.run();
292 396
    checkMatching(graph, mm);
293 397

	
294 398
    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
295 399
    mwm.run();
296 400
    checkWeightedMatching(graph, weight, mwm);
297 401

	
298 402
    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
299 403
    bool perfect = mwpm.run();
300 404

	
301 405
    check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
302 406
          "Perfect matching found");
303 407

	
304 408
    if (perfect) {
0 comments (0 inline)