gravatar
deba@inf.elte.hu
deba@inf.elte.hu
Several improvements in maximum matching algorithms - The interface of MaxMatching is changed to be similar to the weighted algorithms - The internal data structure (the queue implementation and the matching map) is changed in the MaxMatching algorithm, which provides better runtime properties - The Blossom iterators are changed slightly in the weighted matching algorithms - Several documentation improvments - The test files are merged
1 4 0
default
5 files changed with 886 insertions and 1017 deletions:
↑ Collapse diff ↑
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-2008
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
///\brief Maximum matching algorithms in graph.
34
///\brief Maximum matching algorithms in general graphs.
35 35

	
36 36
namespace lemon {
37 37

	
38
  ///\ingroup matching
38
  /// \ingroup matching
39 39
  ///
40
  ///\brief Edmonds' alternating forest maximum matching algorithm.
40
  /// \brief Edmonds' alternating forest maximum matching algorithm.
41 41
  ///
42
  ///This class provides Edmonds' alternating forest matching
43
  ///algorithm. The starting matching (if any) can be passed to the
44
  ///algorithm using some of init functions.
42
  /// This class provides Edmonds' alternating forest matching
43
  /// algorithm. The starting matching (if any) can be passed to the
44
  /// algorithm using some of init functions.
45 45
  ///
46
  ///The dual side of a matching is a map of the nodes to
47
  ///MaxMatching::DecompType, having values \c D, \c A and \c C
48
  ///showing the Gallai-Edmonds decomposition of the digraph. The nodes
49
  ///in \c D induce a digraph with factor-critical components, the nodes
50
  ///in \c A form the barrier, and the nodes in \c C induce a digraph
51
  ///having a perfect matching. This decomposition can be attained by
52
  ///calling \c decomposition() after running the algorithm.
46
  /// The dual side of a matching 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 fractor critical components
53
  /// minus the number of barrier nodes is a lower bound on the
54
  /// unmatched nodes, and if the matching is optimal this bound is
55
  /// tight. This decomposition can be attained by calling \c
56
  /// decomposition() after running the algorithm.
53 57
  ///
54
  ///\param Digraph The graph type the algorithm runs on.
55
  template <typename Graph>
58
  /// \param _Graph The graph type the algorithm runs on.
59
  template <typename _Graph>
56 60
  class MaxMatching {
57

	
58
  protected:
61
  public:
62

	
63
    typedef _Graph Graph;
64
    typedef typename Graph::template NodeMap<typename Graph::Arc>
65
    MatchingMap;
66

	
67
    ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
68
    ///
69
    ///Indicates the Gallai-Edmonds decomposition of the graph, which
70
    ///shows an upper bound on the size of a maximum matching. The
71
    ///nodes with Status \c EVEN/D induce a graph with factor-critical
72
    ///components, the nodes in \c ODD/A form the canonical barrier,
73
    ///and the nodes in \c MATCHED/C induce a graph having a perfect
74
    ///matching.
75
    enum Status {
76
      EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
77
    };
78

	
79
    typedef typename Graph::template NodeMap<Status> StatusMap;
80

	
81
  private:
59 82

	
60 83
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
61 84

	
62
    typedef typename Graph::template NodeMap<int> UFECrossRef;
63
    typedef UnionFindEnum<UFECrossRef> UFE;
64
    typedef std::vector<Node> NV;
65

	
66
    typedef typename Graph::template NodeMap<int> EFECrossRef;
67
    typedef ExtendFindEnum<EFECrossRef> EFE;
85
    typedef UnionFindEnum<IntNodeMap> BlossomSet;
86
    typedef ExtendFindEnum<IntNodeMap> TreeSet;
87
    typedef RangeMap<Node> NodeIntMap;
88
    typedef MatchingMap EarMap;
89
    typedef std::vector<Node> NodeQueue;
90

	
91
    const Graph& _graph;
92
    MatchingMap* _matching;
93
    StatusMap* _status;
94

	
95
    EarMap* _ear;
96

	
97
    IntNodeMap* _blossom_set_index;
98
    BlossomSet* _blossom_set;
99
    NodeIntMap* _blossom_rep;
100

	
101
    IntNodeMap* _tree_set_index;
102
    TreeSet* _tree_set;
103

	
104
    NodeQueue _node_queue;
105
    int _process, _postpone, _last;
106

	
107
    int _node_num;
108

	
109
  private:
110

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

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

	
159
    void processDense(const Node& n) {
160
      _process = _postpone = _last = 0;
161
      _node_queue[_last++] = n;
162

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

	
176
      while (_postpone != _last) {
177
        Node u = _node_queue[_postpone++];
178

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

	
182
          if ((*_status)[v] == EVEN) {
183
            if (_blossom_set->find(u) != _blossom_set->find(v)) {
184
              shrinkOnEdge(a);
185
            }
186
          }
187

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

	
204
    void processSparse(const Node& n) {
205
      _process = _last = 0;
206
      _node_queue[_last++] = n;
207
      while (_process != _last) {
208
        Node u = _node_queue[_process++];
209
        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
210
          Node v = _graph.target(a);
211

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

	
226
    void shrinkOnEdge(const Edge& e) {
227
      Node nca = INVALID;
228

	
229
      {
230
        std::set<Node> left_set, right_set;
231

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

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

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

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

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

	
279
      {
280

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

	
285
        while (base != nca) {
286
          _ear->set(node, arc);
287

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

	
308
      _blossom_rep->set(_blossom_set->find(nca), nca);
309

	
310
      {
311

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

	
316
        while (base != nca) {
317
          _ear->set(node, arc);
318

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

	
339
      _blossom_rep->set(_blossom_set->find(nca), nca);
340
    }
341

	
342

	
343

	
344
    void extendOnArc(const Arc& a) {
345
      Node base = _graph.source(a);
346
      Node odd = _graph.target(a);
347

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

	
358
    }
359

	
360
    void augmentOnArc(const Arc& a) {
361
      Node even = _graph.source(a);
362
      Node odd = _graph.target(a);
363

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

	
366
      _matching->set(odd, _graph.oppositeArc(a));
367
      _status->set(odd, MATCHED);
368

	
369
      Arc arc = (*_matching)[even];
370
      _matching->set(even, a);
371

	
372
      while (arc != INVALID) {
373
        odd = _graph.target(arc);
374
        arc = (*_ear)[odd];
375
        even = _graph.target(arc);
376
        _matching->set(odd, arc);
377
        arc = (*_matching)[even];
378
        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
379
      }
380

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

	
396
    }
68 397

	
69 398
  public:
70 399

	
71
    ///\brief Indicates the Gallai-Edmonds decomposition of the digraph.
400
    /// \brief Constructor
72 401
    ///
73
    ///Indicates the Gallai-Edmonds decomposition of the digraph, which
74
    ///shows an upper bound on the size of a maximum matching. The
75
    ///nodes with DecompType \c D induce a digraph with factor-critical
76
    ///components, the nodes in \c A form the canonical barrier, and the
77
    ///nodes in \c C induce a digraph having a perfect matching.
78
    enum DecompType {
79
      D=0,
80
      A=1,
81
      C=2
82
    };
83

	
84
  protected:
85

	
86
    static const int HEUR_density=2;
87
    const Graph& g;
88
    typename Graph::template NodeMap<Node> _mate;
89
    typename Graph::template NodeMap<DecompType> position;
90

	
91
  public:
92

	
93
    MaxMatching(const Graph& _g)
94
      : g(_g), _mate(_g), position(_g) {}
95

	
96
    ///\brief Sets the actual matching to the empty matching.
402
    /// Constructor.
403
    MaxMatching(const Graph& graph)
404
      : _graph(graph), _matching(0), _status(0), _ear(0),
405
        _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
406
        _tree_set_index(0), _tree_set(0) {}
407

	
408
    ~MaxMatching() {
409
      destroyStructures();
410
    }
411

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

	
417
    /// If you need more control on the execution, first you must call
418
    /// \ref init(), \ref greedyInit() or \ref matchingInit()
419
    /// functions, then you can start the algorithm with the \ref
420
    /// startParse() or startDense() functions.
421

	
422
    ///@{
423

	
424
    /// \brief Sets the actual matching to the empty matching.
97 425
    ///
98
    ///Sets the actual matching to the empty matching.
426
    /// Sets the actual matching to the empty matching.
99 427
    ///
100 428
    void init() {
101
      for(NodeIt v(g); v!=INVALID; ++v) {
102
        _mate.set(v,INVALID);
103
        position.set(v,C);
429
      createStructures();
430
      for(NodeIt n(_graph); n != INVALID; ++n) {
431
        _matching->set(n, INVALID);
432
        _status->set(n, UNMATCHED);
104 433
      }
105 434
    }
106 435

	
107 436
    ///\brief Finds a greedy matching for initial matching.
108 437
    ///
109 438
    ///For initial matchig it finds a maximal greedy matching.
110 439
    void greedyInit() {
111
      for(NodeIt v(g); v!=INVALID; ++v) {
112
        _mate.set(v,INVALID);
113
        position.set(v,C);
440
      createStructures();
441
      for (NodeIt n(_graph); n != INVALID; ++n) {
442
        _matching->set(n, INVALID);
443
        _status->set(n, UNMATCHED);
114 444
      }
115
      for(NodeIt v(g); v!=INVALID; ++v)
116
        if ( _mate[v]==INVALID ) {
117
          for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
118
            Node y=g.runningNode(e);
119
            if ( _mate[y]==INVALID && y!=v ) {
120
              _mate.set(v,y);
121
              _mate.set(y,v);
445
      for (NodeIt n(_graph); n != INVALID; ++n) {
446
        if ((*_matching)[n] == INVALID) {
447
          for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
448
            Node v = _graph.target(a);
449
            if ((*_matching)[v] == INVALID && v != n) {
450
              _matching->set(n, a);
451
              _status->set(n, MATCHED);
452
              _matching->set(v, _graph.oppositeArc(a));
453
              _status->set(v, MATCHED);
122 454
              break;
123 455
            }
124 456
          }
125 457
        }
126
    }
127

	
128
    ///\brief Initialize the matching from each nodes' mate.
129
    ///
130
    ///Initialize the matching from a \c Node valued \c Node map. This
131
    ///map must be \e symmetric, i.e. if \c map[u]==v then \c
132
    ///map[v]==u must hold, and \c uv will be an arc of the initial
133
    ///matching.
134
    template <typename MateMap>
135
    void mateMapInit(MateMap& map) {
136
      for(NodeIt v(g); v!=INVALID; ++v) {
137
        _mate.set(v,map[v]);
138
        position.set(v,C);
139 458
      }
140 459
    }
141 460

	
142
    ///\brief Initialize the matching from a node map with the
143
    ///incident matching arcs.
461

	
462
    /// \brief Initialize the matching from the map containing a matching.
144 463
    ///
145
    ///Initialize the matching from an \c Edge valued \c Node map. \c
146
    ///map[v] must be an \c Edge incident to \c v. This map must have
147
    ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c
148
    ///g.oppositeNode(v,map[v])==u holds, and now some arc joining \c
149
    ///u to \c v will be an arc of the matching.
150
    template<typename MatchingMap>
151
    void matchingMapInit(MatchingMap& map) {
152
      for(NodeIt v(g); v!=INVALID; ++v) {
153
        position.set(v,C);
154
        Edge e=map[v];
155
        if ( e!=INVALID )
156
          _mate.set(v,g.oppositeNode(v,e));
157
        else
158
          _mate.set(v,INVALID);
464
    /// Initialize the matching from a \c bool valued \c Edge map. This
465
    /// map must have the property that there are no two incident edges
466
    /// with true value, ie. it contains a matching.
467
    /// \return %True if the map contains a matching.
468
    template <typename MatchingMap>
469
    bool matchingInit(const MatchingMap& matching) {
470
      createStructures();
471

	
472
      for (NodeIt n(_graph); n != INVALID; ++n) {
473
        _matching->set(n, INVALID);
474
        _status->set(n, UNMATCHED);
159 475
      }
476
      for(EdgeIt e(_graph); e!=INVALID; ++e) {
477
        if (matching[e]) {
478

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

	
484
          Node v = _graph.v(e);
485
          if ((*_matching)[v] != INVALID) return false;
486
          _matching->set(v, _graph.direct(e, false));
487
          _status->set(v, MATCHED);
488
        }
489
      }
490
      return true;
160 491
    }
161 492

	
162
    ///\brief Initialize the matching from the map containing the
163
    ///undirected matching arcs.
493
    /// \brief Starts Edmonds' algorithm
164 494
    ///
165
    ///Initialize the matching from a \c bool valued \c Edge map. This
166
    ///map must have the property that there are no two incident arcs
167
    ///\c e, \c f with \c map[e]==map[f]==true. The arcs \c e with \c
168
    ///map[e]==true form the matching.
169
    template <typename MatchingMap>
170
    void matchingInit(MatchingMap& map) {
171
      for(NodeIt v(g); v!=INVALID; ++v) {
172
        _mate.set(v,INVALID);
173
        position.set(v,C);
174
      }
175
      for(EdgeIt e(g); e!=INVALID; ++e) {
176
        if ( map[e] ) {
177
          Node u=g.u(e);
178
          Node v=g.v(e);
179
          _mate.set(u,v);
180
          _mate.set(v,u);
495
    /// If runs the original Edmonds' algorithm.
496
    void startSparse() {
497
      for(NodeIt n(_graph); n != INVALID; ++n) {
498
        if ((*_status)[n] == UNMATCHED) {
499
          (*_blossom_rep)[_blossom_set->insert(n)] = n;
500
          _tree_set->insert(n);
501
          _status->set(n, EVEN);
502
          processSparse(n);
181 503
        }
182 504
      }
183 505
    }
184 506

	
185

	
186
    ///\brief Runs Edmonds' algorithm.
507
    /// \brief Starts Edmonds' algorithm.
187 508
    ///
188
    ///Runs Edmonds' algorithm for sparse digraphs (number of arcs <
189
    ///2*number of nodes), and a heuristical Edmonds' algorithm with a
190
    ///heuristic of postponing shrinks for dense digraphs.
509
    /// It runs Edmonds' algorithm with a heuristic of postponing
510
    /// shrinks, giving a faster algorithm for dense graphs.
511
    void startDense() {
512
      for(NodeIt n(_graph); n != INVALID; ++n) {
513
        if ((*_status)[n] == UNMATCHED) {
514
          (*_blossom_rep)[_blossom_set->insert(n)] = n;
515
          _tree_set->insert(n);
516
          _status->set(n, EVEN);
517
          processDense(n);
518
        }
519
      }
520
    }
521

	
522

	
523
    /// \brief Runs Edmonds' algorithm
524
    ///
525
    /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
526
    /// or Edmonds' algorithm with a heuristic of
527
    /// postponing shrinks for dense graphs.
191 528
    void run() {
192
      if (countEdges(g) < HEUR_density * countNodes(g)) {
529
      if (countEdges(_graph) < 2 * countNodes(_graph)) {
193 530
        greedyInit();
194 531
        startSparse();
195 532
      } else {
196 533
        init();
197 534
        startDense();
198 535
      }
199 536
    }
200 537

	
201

	
202
    ///\brief Starts Edmonds' algorithm.
203
    ///
204
    ///If runs the original Edmonds' algorithm.
205
    void startSparse() {
206

	
207
      typename Graph::template NodeMap<Node> ear(g,INVALID);
208
      //undefined for the base nodes of the blossoms (i.e. for the
209
      //representative elements of UFE blossom) and for the nodes in C
210

	
211
      UFECrossRef blossom_base(g);
212
      UFE blossom(blossom_base);
213
      NV rep(countNodes(g));
214

	
215
      EFECrossRef tree_base(g);
216
      EFE tree(tree_base);
217

	
218
      //If these UFE's would be members of the class then also
219
      //blossom_base and tree_base should be a member.
220

	
221
      //We build only one tree and the other vertices uncovered by the
222
      //matching belong to C. (They can be considered as singleton
223
      //trees.) If this tree can be augmented or no more
224
      //grow/augmentation/shrink is possible then we return to this
225
      //"for" cycle.
226
      for(NodeIt v(g); v!=INVALID; ++v) {
227
        if (position[v]==C && _mate[v]==INVALID) {
228
          rep[blossom.insert(v)] = v;
229
          tree.insert(v);
230
          position.set(v,D);
231
          normShrink(v, ear, blossom, rep, tree);
232
        }
233
      }
234
    }
235

	
236
    ///\brief Starts Edmonds' algorithm.
237
    ///
238
    ///It runs Edmonds' algorithm with a heuristic of postponing
239
    ///shrinks, giving a faster algorithm for dense digraphs.
240
    void startDense() {
241

	
242
      typename Graph::template NodeMap<Node> ear(g,INVALID);
243
      //undefined for the base nodes of the blossoms (i.e. for the
244
      //representative elements of UFE blossom) and for the nodes in C
245

	
246
      UFECrossRef blossom_base(g);
247
      UFE blossom(blossom_base);
248
      NV rep(countNodes(g));
249

	
250
      EFECrossRef tree_base(g);
251
      EFE tree(tree_base);
252

	
253
      //If these UFE's would be members of the class then also
254
      //blossom_base and tree_base should be a member.
255

	
256
      //We build only one tree and the other vertices uncovered by the
257
      //matching belong to C. (They can be considered as singleton
258
      //trees.) If this tree can be augmented or no more
259
      //grow/augmentation/shrink is possible then we return to this
260
      //"for" cycle.
261
      for(NodeIt v(g); v!=INVALID; ++v) {
262
        if ( position[v]==C && _mate[v]==INVALID ) {
263
          rep[blossom.insert(v)] = v;
264
          tree.insert(v);
265
          position.set(v,D);
266
          lateShrink(v, ear, blossom, rep, tree);
267
        }
268
      }
269
    }
270

	
271

	
538
    /// @}
539

	
540
    /// \name Primal solution
541
    /// Functions for get the primal solution, ie. the matching.
542

	
543
    /// @{
272 544

	
273 545
    ///\brief Returns the size of the actual matching stored.
274 546
    ///
275 547
    ///Returns the size of the actual matching stored. After \ref
276
    ///run() it returns the size of a maximum matching in the digraph.
277
    int size() const {
278
      int s=0;
279
      for(NodeIt v(g); v!=INVALID; ++v) {
280
        if ( _mate[v]!=INVALID ) {
281
          ++s;
548
    ///run() it returns the size of the maximum matching in the graph.
549
    int matchingSize() const {
550
      int size = 0;
551
      for (NodeIt n(_graph); n != INVALID; ++n) {
552
        if ((*_matching)[n] != INVALID) {
553
          ++size;
282 554
        }
283 555
      }
284
      return s/2;
556
      return size / 2;
285 557
    }
286 558

	
559
    /// \brief Returns true when the edge is in the matching.
560
    ///
561
    /// Returns true when the edge is in the matching.
562
    bool matching(const Edge& edge) const {
563
      return edge == (*_matching)[_graph.u(edge)];
564
    }
565

	
566
    /// \brief Returns the matching edge incident to the given node.
567
    ///
568
    /// Returns the matching edge of a \c node in the actual matching or
569
    /// INVALID if the \c node is not covered by the actual matching.
570
    Arc matching(const Node& n) const {
571
      return (*_matching)[n];
572
    }
287 573

	
288 574
    ///\brief Returns the mate of a node in the actual matching.
289 575
    ///
290
    ///Returns the mate of a \c node in the actual matching.
291
    ///Returns INVALID if the \c node is not covered by the actual matching.
292
    Node mate(const Node& node) const {
293
      return _mate[node];
576
    ///Returns the mate of a \c node in the actual matching or
577
    ///INVALID if the \c node is not covered by the actual matching.
578
    Node mate(const Node& n) const {
579
      return (*_matching)[n] != INVALID ?
580
        _graph.target((*_matching)[n]) : INVALID;
294 581
    }
295 582

	
296
    ///\brief Returns the matching arc incident to the given node.
297
    ///
298
    ///Returns the matching arc of a \c node in the actual matching.
299
    ///Returns INVALID if the \c node is not covered by the actual matching.
300
    Edge matchingArc(const Node& node) const {
301
      if (_mate[node] == INVALID) return INVALID;
302
      Node n = node < _mate[node] ? node : _mate[node];
303
      for (IncEdgeIt e(g, n); e != INVALID; ++e) {
304
        if (g.oppositeNode(n, e) == _mate[n]) {
305
          return e;
306
        }
307
      }
308
      return INVALID;
309
    }
583
    /// @}
584

	
585
    /// \name Dual solution
586
    /// Functions for get the dual solution, ie. the decomposition.
587

	
588
    /// @{
310 589

	
311 590
    /// \brief Returns the class of the node in the Edmonds-Gallai
312 591
    /// decomposition.
313 592
    ///
314 593
    /// Returns the class of the node in the Edmonds-Gallai
315 594
    /// decomposition.
316
    DecompType decomposition(const Node& n) {
317
      return position[n] == A;
595
    Status decomposition(const Node& n) const {
596
      return (*_status)[n];
318 597
    }
319 598

	
320 599
    /// \brief Returns true when the node is in the barrier.
321 600
    ///
322 601
    /// Returns true when the node is in the barrier.
323
    bool barrier(const Node& n) {
324
      return position[n] == A;
602
    bool barrier(const Node& n) const {
603
      return (*_status)[n] == ODD;
325 604
    }
326 605

	
327
    ///\brief Gives back the matching in a \c Node of mates.
328
    ///
329
    ///Writes the stored matching to a \c Node valued \c Node map. The
330
    ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
331
    ///map[v]==u will hold, and now \c uv is an arc of the matching.
332
    template <typename MateMap>
333
    void mateMap(MateMap& map) const {
334
      for(NodeIt v(g); v!=INVALID; ++v) {
335
        map.set(v,_mate[v]);
336
      }
337
    }
338

	
339
    ///\brief Gives back the matching in an \c Edge valued \c Node
340
    ///map.
341
    ///
342
    ///Writes the stored matching to an \c Edge valued \c Node
343
    ///map. \c map[v] will be an \c Edge incident to \c v. This
344
    ///map will have the property that if \c g.oppositeNode(u,map[u])
345
    ///== v then \c map[u]==map[v] holds, and now this arc is an arc
346
    ///of the matching.
347
    template <typename MatchingMap>
348
    void matchingMap(MatchingMap& map)  const {
349
      typename Graph::template NodeMap<bool> todo(g,true);
350
      for(NodeIt v(g); v!=INVALID; ++v) {
351
        if (_mate[v]!=INVALID && v < _mate[v]) {
352
          Node u=_mate[v];
353
          for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
354
            if ( g.runningNode(e) == u ) {
355
              map.set(u,e);
356
              map.set(v,e);
357
              todo.set(u,false);
358
              todo.set(v,false);
359
              break;
360
            }
361
          }
362
        }
363
      }
364
    }
365

	
366

	
367
    ///\brief Gives back the matching in a \c bool valued \c Edge
368
    ///map.
369
    ///
370
    ///Writes the matching stored to a \c bool valued \c Arc
371
    ///map. This map will have the property that there are no two
372
    ///incident arcs \c e, \c f with \c map[e]==map[f]==true. The
373
    ///arcs \c e with \c map[e]==true form the matching.
374
    template<typename MatchingMap>
375
    void matching(MatchingMap& map) const {
376
      for(EdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
377

	
378
      typename Graph::template NodeMap<bool> todo(g,true);
379
      for(NodeIt v(g); v!=INVALID; ++v) {
380
        if ( todo[v] && _mate[v]!=INVALID ) {
381
          Node u=_mate[v];
382
          for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
383
            if ( g.runningNode(e) == u ) {
384
              map.set(e,true);
385
              todo.set(u,false);
386
              todo.set(v,false);
387
              break;
388
            }
389
          }
390
        }
391
      }
392
    }
393

	
394

	
395
    ///\brief Returns the canonical decomposition of the digraph after running
396
    ///the algorithm.
397
    ///
398
    ///After calling any run methods of the class, it writes the
399
    ///Gallai-Edmonds canonical decomposition of the digraph. \c map
400
    ///must be a node map of \ref DecompType 's.
401
    template <typename DecompositionMap>
402
    void decomposition(DecompositionMap& map) const {
403
      for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]);
404
    }
405

	
406
    ///\brief Returns a barrier on the nodes.
407
    ///
408
    ///After calling any run methods of the class, it writes a
409
    ///canonical barrier on the nodes. The odd component number of the
410
    ///remaining digraph minus the barrier size is a lower bound for the
411
    ///uncovered nodes in the digraph. The \c map must be a node map of
412
    ///bools.
413
    template <typename BarrierMap>
414
    void barrier(BarrierMap& barrier) {
415
      for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A);
416
    }
417

	
418
  private:
419

	
420

	
421
    void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,
422
                    UFE& blossom, NV& rep, EFE& tree) {
423
      //We have one tree which we grow, and also shrink but only if it
424
      //cannot be postponed. If we augment then we return to the "for"
425
      //cycle of runEdmonds().
426

	
427
      std::queue<Node> Q;   //queue of the totally unscanned nodes
428
      Q.push(v);
429
      std::queue<Node> R;
430
      //queue of the nodes which must be scanned for a possible shrink
431

	
432
      while ( !Q.empty() ) {
433
        Node x=Q.front();
434
        Q.pop();
435
        for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
436
          Node y=g.runningNode(e);
437
          //growOrAugment grows if y is covered by the matching and
438
          //augments if not. In this latter case it returns 1.
439
          if (position[y]==C &&
440
              growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
441
        }
442
        R.push(x);
443
      }
444

	
445
      while ( !R.empty() ) {
446
        Node x=R.front();
447
        R.pop();
448

	
449
        for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
450
          Node y=g.runningNode(e);
451

	
452
          if ( position[y] == D && blossom.find(x) != blossom.find(y) )
453
            //Recall that we have only one tree.
454
            shrink( x, y, ear, blossom, rep, tree, Q);
455

	
456
          while ( !Q.empty() ) {
457
            Node z=Q.front();
458
            Q.pop();
459
            for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
460
              Node w=g.runningNode(f);
461
              //growOrAugment grows if y is covered by the matching and
462
              //augments if not. In this latter case it returns 1.
463
              if (position[w]==C &&
464
                  growOrAugment(w, z, ear, blossom, rep, tree, Q)) return;
465
            }
466
            R.push(z);
467
          }
468
        } //for e
469
      } // while ( !R.empty() )
470
    }
471

	
472
    void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,
473
                    UFE& blossom, NV& rep, EFE& tree) {
474
      //We have one tree, which we grow and shrink. If we augment then we
475
      //return to the "for" cycle of runEdmonds().
476

	
477
      std::queue<Node> Q;   //queue of the unscanned nodes
478
      Q.push(v);
479
      while ( !Q.empty() ) {
480

	
481
        Node x=Q.front();
482
        Q.pop();
483

	
484
        for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
485
          Node y=g.runningNode(e);
486

	
487
          switch ( position[y] ) {
488
          case D:          //x and y must be in the same tree
489
            if ( blossom.find(x) != blossom.find(y))
490
              //x and y are in the same tree
491
              shrink(x, y, ear, blossom, rep, tree, Q);
492
            break;
493
          case C:
494
            //growOrAugment grows if y is covered by the matching and
495
            //augments if not. In this latter case it returns 1.
496
            if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
497
            break;
498
          default: break;
499
          }
500
        }
501
      }
502
    }
503

	
504
    void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear,
505
                UFE& blossom, NV& rep, EFE& tree,std::queue<Node>& Q) {
506
      //x and y are the two adjacent vertices in two blossoms.
507

	
508
      typename Graph::template NodeMap<bool> path(g,false);
509

	
510
      Node b=rep[blossom.find(x)];
511
      path.set(b,true);
512
      b=_mate[b];
513
      while ( b!=INVALID ) {
514
        b=rep[blossom.find(ear[b])];
515
        path.set(b,true);
516
        b=_mate[b];
517
      } //we go until the root through bases of blossoms and odd vertices
518

	
519
      Node top=y;
520
      Node middle=rep[blossom.find(top)];
521
      Node bottom=x;
522
      while ( !path[middle] )
523
        shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
524
      //Until we arrive to a node on the path, we update blossom, tree
525
      //and the positions of the odd nodes.
526

	
527
      Node base=middle;
528
      top=x;
529
      middle=rep[blossom.find(top)];
530
      bottom=y;
531
      Node blossom_base=rep[blossom.find(base)];
532
      while ( middle!=blossom_base )
533
        shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
534
      //Until we arrive to a node on the path, we update blossom, tree
535
      //and the positions of the odd nodes.
536

	
537
      rep[blossom.find(base)] = base;
538
    }
539

	
540
    void shrinkStep(Node& top, Node& middle, Node& bottom,
541
                    typename Graph::template NodeMap<Node>& ear,
542
                    UFE& blossom, NV& rep, EFE& tree, std::queue<Node>& Q) {
543
      //We traverse a blossom and update everything.
544

	
545
      ear.set(top,bottom);
546
      Node t=top;
547
      while ( t!=middle ) {
548
        Node u=_mate[t];
549
        t=ear[u];
550
        ear.set(t,u);
551
      }
552
      bottom=_mate[middle];
553
      position.set(bottom,D);
554
      Q.push(bottom);
555
      top=ear[bottom];
556
      Node oldmiddle=middle;
557
      middle=rep[blossom.find(top)];
558
      tree.erase(bottom);
559
      tree.erase(oldmiddle);
560
      blossom.insert(bottom);
561
      blossom.join(bottom, oldmiddle);
562
      blossom.join(top, oldmiddle);
563
    }
564

	
565

	
566

	
567
    bool growOrAugment(Node& y, Node& x, typename Graph::template
568
                       NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree,
569
                       std::queue<Node>& Q) {
570
      //x is in a blossom in the tree, y is outside. If y is covered by
571
      //the matching we grow, otherwise we augment. In this case we
572
      //return 1.
573

	
574
      if ( _mate[y]!=INVALID ) {       //grow
575
        ear.set(y,x);
576
        Node w=_mate[y];
577
        rep[blossom.insert(w)] = w;
578
        position.set(y,A);
579
        position.set(w,D);
580
        int t = tree.find(rep[blossom.find(x)]);
581
        tree.insert(y,t);
582
        tree.insert(w,t);
583
        Q.push(w);
584
      } else {                      //augment
585
        augment(x, ear, blossom, rep, tree);
586
        _mate.set(x,y);
587
        _mate.set(y,x);
588
        return true;
589
      }
590
      return false;
591
    }
592

	
593
    void augment(Node x, typename Graph::template NodeMap<Node>& ear,
594
                 UFE& blossom, NV& rep, EFE& tree) {
595
      Node v=_mate[x];
596
      while ( v!=INVALID ) {
597

	
598
        Node u=ear[v];
599
        _mate.set(v,u);
600
        Node tmp=v;
601
        v=_mate[u];
602
        _mate.set(u,tmp);
603
      }
604
      int y = tree.find(rep[blossom.find(x)]);
605
      for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {
606
        if ( position[tit] == D ) {
607
          int b = blossom.find(tit);
608
          for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) {
609
            position.set(bit, C);
610
          }
611
          blossom.eraseClass(b);
612
        } else position.set(tit, C);
613
      }
614
      tree.eraseClass(y);
615

	
616
    }
606
    /// @}
617 607

	
618 608
  };
619 609

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

	
663 656
    typedef _Graph Graph;
664 657
    typedef _WeightMap WeightMap;
665 658
    typedef typename WeightMap::Value Value;
666 659

	
667 660
    /// \brief Scaling factor for dual solution
668 661
    ///
669 662
    /// Scaling factor for dual solution, it is equal to 4 or 1
670 663
    /// according to the value type.
671 664
    static const int dualScale =
672 665
      std::numeric_limits<Value>::is_integer ? 4 : 1;
673 666

	
674 667
    typedef typename Graph::template NodeMap<typename Graph::Arc>
675 668
    MatchingMap;
676 669

	
677 670
  private:
678 671

	
679 672
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
680 673

	
681 674
    typedef typename Graph::template NodeMap<Value> NodePotential;
682 675
    typedef std::vector<Node> BlossomNodeList;
683 676

	
684 677
    struct BlossomVariable {
685 678
      int begin, end;
686 679
      Value value;
687 680

	
688 681
      BlossomVariable(int _begin, int _end, Value _value)
689 682
        : begin(_begin), end(_end), value(_value) {}
690 683

	
691 684
    };
692 685

	
693 686
    typedef std::vector<BlossomVariable> BlossomPotential;
694 687

	
695 688
    const Graph& _graph;
696 689
    const WeightMap& _weight;
697 690

	
698 691
    MatchingMap* _matching;
699 692

	
700 693
    NodePotential* _node_potential;
701 694

	
702 695
    BlossomPotential _blossom_potential;
703 696
    BlossomNodeList _blossom_node_list;
704 697

	
705 698
    int _node_num;
706 699
    int _blossom_num;
707 700

	
708
    typedef typename Graph::template NodeMap<int> NodeIntMap;
709
    typedef typename Graph::template ArcMap<int> ArcIntMap;
710
    typedef typename Graph::template EdgeMap<int> EdgeIntMap;
711 701
    typedef RangeMap<int> IntIntMap;
712 702

	
713 703
    enum Status {
714 704
      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
715 705
    };
716 706

	
717
    typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
707
    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
718 708
    struct BlossomData {
719 709
      int tree;
720 710
      Status status;
721 711
      Arc pred, next;
722 712
      Value pot, offset;
723 713
      Node base;
724 714
    };
725 715

	
726
    NodeIntMap *_blossom_index;
716
    IntNodeMap *_blossom_index;
727 717
    BlossomSet *_blossom_set;
728 718
    RangeMap<BlossomData>* _blossom_data;
729 719

	
730
    NodeIntMap *_node_index;
731
    ArcIntMap *_node_heap_index;
720
    IntNodeMap *_node_index;
721
    IntArcMap *_node_heap_index;
732 722

	
733 723
    struct NodeData {
734 724

	
735
      NodeData(ArcIntMap& node_heap_index)
725
      NodeData(IntArcMap& node_heap_index)
736 726
        : heap(node_heap_index) {}
737 727

	
738 728
      int blossom;
739 729
      Value pot;
740
      BinHeap<Value, ArcIntMap> heap;
730
      BinHeap<Value, IntArcMap> heap;
741 731
      std::map<int, Arc> heap_index;
742 732

	
743 733
      int tree;
744 734
    };
745 735

	
746 736
    RangeMap<NodeData>* _node_data;
747 737

	
748 738
    typedef ExtendFindEnum<IntIntMap> TreeSet;
749 739

	
750 740
    IntIntMap *_tree_set_index;
751 741
    TreeSet *_tree_set;
752 742

	
753
    NodeIntMap *_delta1_index;
754
    BinHeap<Value, NodeIntMap> *_delta1;
743
    IntNodeMap *_delta1_index;
744
    BinHeap<Value, IntNodeMap> *_delta1;
755 745

	
756 746
    IntIntMap *_delta2_index;
757 747
    BinHeap<Value, IntIntMap> *_delta2;
758 748

	
759
    EdgeIntMap *_delta3_index;
760
    BinHeap<Value, EdgeIntMap> *_delta3;
749
    IntEdgeMap *_delta3_index;
750
    BinHeap<Value, IntEdgeMap> *_delta3;
761 751

	
762 752
    IntIntMap *_delta4_index;
763 753
    BinHeap<Value, IntIntMap> *_delta4;
764 754

	
765 755
    Value _delta_sum;
766 756

	
767 757
    void createStructures() {
768 758
      _node_num = countNodes(_graph);
769 759
      _blossom_num = _node_num * 3 / 2;
770 760

	
771 761
      if (!_matching) {
772 762
        _matching = new MatchingMap(_graph);
773 763
      }
774 764
      if (!_node_potential) {
775 765
        _node_potential = new NodePotential(_graph);
776 766
      }
777 767
      if (!_blossom_set) {
778
        _blossom_index = new NodeIntMap(_graph);
768
        _blossom_index = new IntNodeMap(_graph);
779 769
        _blossom_set = new BlossomSet(*_blossom_index);
780 770
        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
781 771
      }
782 772

	
783 773
      if (!_node_index) {
784
        _node_index = new NodeIntMap(_graph);
785
        _node_heap_index = new ArcIntMap(_graph);
774
        _node_index = new IntNodeMap(_graph);
775
        _node_heap_index = new IntArcMap(_graph);
786 776
        _node_data = new RangeMap<NodeData>(_node_num,
787 777
                                              NodeData(*_node_heap_index));
788 778
      }
789 779

	
790 780
      if (!_tree_set) {
791 781
        _tree_set_index = new IntIntMap(_blossom_num);
792 782
        _tree_set = new TreeSet(*_tree_set_index);
793 783
      }
794 784
      if (!_delta1) {
795
        _delta1_index = new NodeIntMap(_graph);
796
        _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);
785
        _delta1_index = new IntNodeMap(_graph);
786
        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
797 787
      }
798 788
      if (!_delta2) {
799 789
        _delta2_index = new IntIntMap(_blossom_num);
800 790
        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
801 791
      }
802 792
      if (!_delta3) {
803
        _delta3_index = new EdgeIntMap(_graph);
804
        _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
793
        _delta3_index = new IntEdgeMap(_graph);
794
        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
805 795
      }
806 796
      if (!_delta4) {
807 797
        _delta4_index = new IntIntMap(_blossom_num);
808 798
        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
809 799
      }
810 800
    }
811 801

	
812 802
    void destroyStructures() {
813 803
      _node_num = countNodes(_graph);
814 804
      _blossom_num = _node_num * 3 / 2;
815 805

	
816 806
      if (_matching) {
817 807
        delete _matching;
818 808
      }
819 809
      if (_node_potential) {
820 810
        delete _node_potential;
821 811
      }
822 812
      if (_blossom_set) {
823 813
        delete _blossom_index;
824 814
        delete _blossom_set;
825 815
        delete _blossom_data;
826 816
      }
827 817

	
828 818
      if (_node_index) {
829 819
        delete _node_index;
830 820
        delete _node_heap_index;
831 821
        delete _node_data;
832 822
      }
833 823

	
834 824
      if (_tree_set) {
835 825
        delete _tree_set_index;
836 826
        delete _tree_set;
837 827
      }
838 828
      if (_delta1) {
839 829
        delete _delta1_index;
840 830
        delete _delta1;
841 831
      }
842 832
      if (_delta2) {
843 833
        delete _delta2_index;
844 834
        delete _delta2;
845 835
      }
846 836
      if (_delta3) {
847 837
        delete _delta3_index;
848 838
        delete _delta3;
849 839
      }
850 840
      if (_delta4) {
851 841
        delete _delta4_index;
852 842
        delete _delta4;
... ...
@@ -1221,141 +1211,141 @@
1221 1211
    void alternatePath(int even, int tree) {
1222 1212
      int odd;
1223 1213

	
1224 1214
      evenToMatched(even, tree);
1225 1215
      (*_blossom_data)[even].status = MATCHED;
1226 1216

	
1227 1217
      while ((*_blossom_data)[even].pred != INVALID) {
1228 1218
        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
1229 1219
        (*_blossom_data)[odd].status = MATCHED;
1230 1220
        oddToMatched(odd);
1231 1221
        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1232 1222

	
1233 1223
        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
1234 1224
        (*_blossom_data)[even].status = MATCHED;
1235 1225
        evenToMatched(even, tree);
1236 1226
        (*_blossom_data)[even].next =
1237 1227
          _graph.oppositeArc((*_blossom_data)[odd].pred);
1238 1228
      }
1239 1229

	
1240 1230
    }
1241 1231

	
1242 1232
    void destroyTree(int tree) {
1243 1233
      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1244 1234
        if ((*_blossom_data)[b].status == EVEN) {
1245 1235
          (*_blossom_data)[b].status = MATCHED;
1246 1236
          evenToMatched(b, tree);
1247 1237
        } else if ((*_blossom_data)[b].status == ODD) {
1248 1238
          (*_blossom_data)[b].status = MATCHED;
1249 1239
          oddToMatched(b);
1250 1240
        }
1251 1241
      }
1252 1242
      _tree_set->eraseClass(tree);
1253 1243
    }
1254 1244

	
1255 1245

	
1256 1246
    void unmatchNode(const Node& node) {
1257 1247
      int blossom = _blossom_set->find(node);
1258 1248
      int tree = _tree_set->find(blossom);
1259 1249

	
1260 1250
      alternatePath(blossom, tree);
1261 1251
      destroyTree(tree);
1262 1252

	
1263 1253
      (*_blossom_data)[blossom].status = UNMATCHED;
1264 1254
      (*_blossom_data)[blossom].base = node;
1265 1255
      matchedToUnmatched(blossom);
1266 1256
    }
1267 1257

	
1268 1258

	
1269
    void augmentOnArc(const Edge& arc) {
1270

	
1271
      int left = _blossom_set->find(_graph.u(arc));
1272
      int right = _blossom_set->find(_graph.v(arc));
1259
    void augmentOnEdge(const Edge& edge) {
1260

	
1261
      int left = _blossom_set->find(_graph.u(edge));
1262
      int right = _blossom_set->find(_graph.v(edge));
1273 1263

	
1274 1264
      if ((*_blossom_data)[left].status == EVEN) {
1275 1265
        int left_tree = _tree_set->find(left);
1276 1266
        alternatePath(left, left_tree);
1277 1267
        destroyTree(left_tree);
1278 1268
      } else {
1279 1269
        (*_blossom_data)[left].status = MATCHED;
1280 1270
        unmatchedToMatched(left);
1281 1271
      }
1282 1272

	
1283 1273
      if ((*_blossom_data)[right].status == EVEN) {
1284 1274
        int right_tree = _tree_set->find(right);
1285 1275
        alternatePath(right, right_tree);
1286 1276
        destroyTree(right_tree);
1287 1277
      } else {
1288 1278
        (*_blossom_data)[right].status = MATCHED;
1289 1279
        unmatchedToMatched(right);
1290 1280
      }
1291 1281

	
1292
      (*_blossom_data)[left].next = _graph.direct(arc, true);
1293
      (*_blossom_data)[right].next = _graph.direct(arc, false);
1282
      (*_blossom_data)[left].next = _graph.direct(edge, true);
1283
      (*_blossom_data)[right].next = _graph.direct(edge, false);
1294 1284
    }
1295 1285

	
1296 1286
    void extendOnArc(const Arc& arc) {
1297 1287
      int base = _blossom_set->find(_graph.target(arc));
1298 1288
      int tree = _tree_set->find(base);
1299 1289

	
1300 1290
      int odd = _blossom_set->find(_graph.source(arc));
1301 1291
      _tree_set->insert(odd, tree);
1302 1292
      (*_blossom_data)[odd].status = ODD;
1303 1293
      matchedToOdd(odd);
1304 1294
      (*_blossom_data)[odd].pred = arc;
1305 1295

	
1306 1296
      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
1307 1297
      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1308 1298
      _tree_set->insert(even, tree);
1309 1299
      (*_blossom_data)[even].status = EVEN;
1310 1300
      matchedToEven(even, tree);
1311 1301
    }
1312 1302

	
1313
    void shrinkOnArc(const Edge& edge, int tree) {
1303
    void shrinkOnEdge(const Edge& edge, int tree) {
1314 1304
      int nca = -1;
1315 1305
      std::vector<int> left_path, right_path;
1316 1306

	
1317 1307
      {
1318 1308
        std::set<int> left_set, right_set;
1319 1309
        int left = _blossom_set->find(_graph.u(edge));
1320 1310
        left_path.push_back(left);
1321 1311
        left_set.insert(left);
1322 1312

	
1323 1313
        int right = _blossom_set->find(_graph.v(edge));
1324 1314
        right_path.push_back(right);
1325 1315
        right_set.insert(right);
1326 1316

	
1327 1317
        while (true) {
1328 1318

	
1329 1319
          if ((*_blossom_data)[left].pred == INVALID) break;
1330 1320

	
1331 1321
          left =
1332 1322
            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1333 1323
          left_path.push_back(left);
1334 1324
          left =
1335 1325
            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1336 1326
          left_path.push_back(left);
1337 1327

	
1338 1328
          left_set.insert(left);
1339 1329

	
1340 1330
          if (right_set.find(left) != right_set.end()) {
1341 1331
            nca = left;
1342 1332
            break;
1343 1333
          }
1344 1334

	
1345 1335
          if ((*_blossom_data)[right].pred == INVALID) break;
1346 1336

	
1347 1337
          right =
1348 1338
            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1349 1339
          right_path.push_back(right);
1350 1340
          right =
1351 1341
            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1352 1342
          right_path.push_back(right);
1353 1343

	
1354 1344
          right_set.insert(right);
1355 1345

	
1356 1346
          if (left_set.find(right) != left_set.end()) {
1357 1347
            nca = right;
1358 1348
            break;
1359 1349
          }
1360 1350

	
1361 1351
        }
... ...
@@ -1607,97 +1597,97 @@
1607 1597
          Arc matching = (*_blossom_data)[blossoms[i]].next;
1608 1598
          Node base = _graph.source(matching);
1609 1599
          extractBlossom(blossoms[i], base, matching);
1610 1600
        } else {
1611 1601
          Node base = (*_blossom_data)[blossoms[i]].base;
1612 1602
          extractBlossom(blossoms[i], base, INVALID);
1613 1603
        }
1614 1604
      }
1615 1605
    }
1616 1606

	
1617 1607
  public:
1618 1608

	
1619 1609
    /// \brief Constructor
1620 1610
    ///
1621 1611
    /// Constructor.
1622 1612
    MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1623 1613
      : _graph(graph), _weight(weight), _matching(0),
1624 1614
        _node_potential(0), _blossom_potential(), _blossom_node_list(),
1625 1615
        _node_num(0), _blossom_num(0),
1626 1616

	
1627 1617
        _blossom_index(0), _blossom_set(0), _blossom_data(0),
1628 1618
        _node_index(0), _node_heap_index(0), _node_data(0),
1629 1619
        _tree_set_index(0), _tree_set(0),
1630 1620

	
1631 1621
        _delta1_index(0), _delta1(0),
1632 1622
        _delta2_index(0), _delta2(0),
1633 1623
        _delta3_index(0), _delta3(0),
1634 1624
        _delta4_index(0), _delta4(0),
1635 1625

	
1636 1626
        _delta_sum() {}
1637 1627

	
1638 1628
    ~MaxWeightedMatching() {
1639 1629
      destroyStructures();
1640 1630
    }
1641 1631

	
1642 1632
    /// \name Execution control
1643 1633
    /// The simplest way to execute the algorithm is to use the member
1644 1634
    /// \c run() member function.
1645 1635

	
1646 1636
    ///@{
1647 1637

	
1648 1638
    /// \brief Initialize the algorithm
1649 1639
    ///
1650 1640
    /// Initialize the algorithm
1651 1641
    void init() {
1652 1642
      createStructures();
1653 1643

	
1654 1644
      for (ArcIt e(_graph); e != INVALID; ++e) {
1655
        _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
1645
        _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
1656 1646
      }
1657 1647
      for (NodeIt n(_graph); n != INVALID; ++n) {
1658 1648
        _delta1_index->set(n, _delta1->PRE_HEAP);
1659 1649
      }
1660 1650
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1661 1651
        _delta3_index->set(e, _delta3->PRE_HEAP);
1662 1652
      }
1663 1653
      for (int i = 0; i < _blossom_num; ++i) {
1664 1654
        _delta2_index->set(i, _delta2->PRE_HEAP);
1665 1655
        _delta4_index->set(i, _delta4->PRE_HEAP);
1666 1656
      }
1667 1657

	
1668 1658
      int index = 0;
1669 1659
      for (NodeIt n(_graph); n != INVALID; ++n) {
1670 1660
        Value max = 0;
1671 1661
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1672 1662
          if (_graph.target(e) == n) continue;
1673 1663
          if ((dualScale * _weight[e]) / 2 > max) {
1674 1664
            max = (dualScale * _weight[e]) / 2;
1675 1665
          }
1676 1666
        }
1677 1667
        _node_index->set(n, index);
1678 1668
        (*_node_data)[index].pot = max;
1679 1669
        _delta1->push(n, max);
1680 1670
        int blossom =
1681 1671
          _blossom_set->insert(n, std::numeric_limits<Value>::max());
1682 1672

	
1683 1673
        _tree_set->insert(blossom);
1684 1674

	
1685 1675
        (*_blossom_data)[blossom].status = EVEN;
1686 1676
        (*_blossom_data)[blossom].pred = INVALID;
1687 1677
        (*_blossom_data)[blossom].next = INVALID;
1688 1678
        (*_blossom_data)[blossom].pot = 0;
1689 1679
        (*_blossom_data)[blossom].offset = 0;
1690 1680
        ++index;
1691 1681
      }
1692 1682
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1693 1683
        int si = (*_node_index)[_graph.u(e)];
1694 1684
        int ti = (*_node_index)[_graph.v(e)];
1695 1685
        if (_graph.u(e) != _graph.v(e)) {
1696 1686
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1697 1687
                            dualScale * _weight[e]) / 2);
1698 1688
        }
1699 1689
      }
1700 1690
    }
1701 1691

	
1702 1692
    /// \brief Starts the algorithm
1703 1693
    ///
... ...
@@ -1724,456 +1714,464 @@
1724 1714
        _delta_sum = d1; OpType ot = D1;
1725 1715
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1726 1716
        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1727 1717
        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1728 1718

	
1729 1719

	
1730 1720
        switch (ot) {
1731 1721
        case D1:
1732 1722
          {
1733 1723
            Node n = _delta1->top();
1734 1724
            unmatchNode(n);
1735 1725
            --unmatched;
1736 1726
          }
1737 1727
          break;
1738 1728
        case D2:
1739 1729
          {
1740 1730
            int blossom = _delta2->top();
1741 1731
            Node n = _blossom_set->classTop(blossom);
1742 1732
            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
1743 1733
            extendOnArc(e);
1744 1734
          }
1745 1735
          break;
1746 1736
        case D3:
1747 1737
          {
1748 1738
            Edge e = _delta3->top();
1749 1739

	
1750 1740
            int left_blossom = _blossom_set->find(_graph.u(e));
1751 1741
            int right_blossom = _blossom_set->find(_graph.v(e));
1752 1742

	
1753 1743
            if (left_blossom == right_blossom) {
1754 1744
              _delta3->pop();
1755 1745
            } else {
1756 1746
              int left_tree;
1757 1747
              if ((*_blossom_data)[left_blossom].status == EVEN) {
1758 1748
                left_tree = _tree_set->find(left_blossom);
1759 1749
              } else {
1760 1750
                left_tree = -1;
1761 1751
                ++unmatched;
1762 1752
              }
1763 1753
              int right_tree;
1764 1754
              if ((*_blossom_data)[right_blossom].status == EVEN) {
1765 1755
                right_tree = _tree_set->find(right_blossom);
1766 1756
              } else {
1767 1757
                right_tree = -1;
1768 1758
                ++unmatched;
1769 1759
              }
1770 1760

	
1771 1761
              if (left_tree == right_tree) {
1772
                shrinkOnArc(e, left_tree);
1762
                shrinkOnEdge(e, left_tree);
1773 1763
              } else {
1774
                augmentOnArc(e);
1764
                augmentOnEdge(e);
1775 1765
                unmatched -= 2;
1776 1766
              }
1777 1767
            }
1778 1768
          } break;
1779 1769
        case D4:
1780 1770
          splitBlossom(_delta4->top());
1781 1771
          break;
1782 1772
        }
1783 1773
      }
1784 1774
      extractMatching();
1785 1775
    }
1786 1776

	
1787 1777
    /// \brief Runs %MaxWeightedMatching algorithm.
1788 1778
    ///
1789 1779
    /// This method runs the %MaxWeightedMatching algorithm.
1790 1780
    ///
1791 1781
    /// \note mwm.run() is just a shortcut of the following code.
1792 1782
    /// \code
1793 1783
    ///   mwm.init();
1794 1784
    ///   mwm.start();
1795 1785
    /// \endcode
1796 1786
    void run() {
1797 1787
      init();
1798 1788
      start();
1799 1789
    }
1800 1790

	
1801 1791
    /// @}
1802 1792

	
1803 1793
    /// \name Primal solution
1804 1794
    /// Functions for get the primal solution, ie. the matching.
1805 1795

	
1806 1796
    /// @{
1807 1797

	
1808 1798
    /// \brief Returns the matching value.
1809 1799
    ///
1810 1800
    /// Returns the matching value.
1811 1801
    Value matchingValue() const {
1812 1802
      Value sum = 0;
1813 1803
      for (NodeIt n(_graph); n != INVALID; ++n) {
1814 1804
        if ((*_matching)[n] != INVALID) {
1815 1805
          sum += _weight[(*_matching)[n]];
1816 1806
        }
1817 1807
      }
1818 1808
      return sum /= 2;
1819 1809
    }
1820 1810

	
1821
    /// \brief Returns true when the arc is in the matching.
1811
    /// \brief Returns the cardinality of the matching.
1822 1812
    ///
1823
    /// Returns true when the arc is in the matching.
1824
    bool matching(const Edge& arc) const {
1825
      return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
1813
    /// Returns the cardinality of the matching.
1814
    int matchingSize() const {
1815
      int num = 0;
1816
      for (NodeIt n(_graph); n != INVALID; ++n) {
1817
        if ((*_matching)[n] != INVALID) {
1818
          ++num;
1819
        }
1820
      }
1821
      return num /= 2;
1822
    }
1823

	
1824
    /// \brief Returns true when the edge is in the matching.
1825
    ///
1826
    /// Returns true when the edge is in the matching.
1827
    bool matching(const Edge& edge) const {
1828
      return edge == (*_matching)[_graph.u(edge)];
1826 1829
    }
1827 1830

	
1828 1831
    /// \brief Returns the incident matching arc.
1829 1832
    ///
1830 1833
    /// Returns the incident matching arc from given node. If the
1831 1834
    /// node is not matched then it gives back \c INVALID.
1832 1835
    Arc matching(const Node& node) const {
1833 1836
      return (*_matching)[node];
1834 1837
    }
1835 1838

	
1836 1839
    /// \brief Returns the mate of the node.
1837 1840
    ///
1838 1841
    /// Returns the adjancent node in a mathcing arc. If the node is
1839 1842
    /// not matched then it gives back \c INVALID.
1840 1843
    Node mate(const Node& node) const {
1841 1844
      return (*_matching)[node] != INVALID ?
1842 1845
        _graph.target((*_matching)[node]) : INVALID;
1843 1846
    }
1844 1847

	
1845 1848
    /// @}
1846 1849

	
1847 1850
    /// \name Dual solution
1848 1851
    /// Functions for get the dual solution.
1849 1852

	
1850 1853
    /// @{
1851 1854

	
1852 1855
    /// \brief Returns the value of the dual solution.
1853 1856
    ///
1854 1857
    /// Returns the value of the dual solution. It should be equal to
1855 1858
    /// the primal value scaled by \ref dualScale "dual scale".
1856 1859
    Value dualValue() const {
1857 1860
      Value sum = 0;
1858 1861
      for (NodeIt n(_graph); n != INVALID; ++n) {
1859 1862
        sum += nodeValue(n);
1860 1863
      }
1861 1864
      for (int i = 0; i < blossomNum(); ++i) {
1862 1865
        sum += blossomValue(i) * (blossomSize(i) / 2);
1863 1866
      }
1864 1867
      return sum;
1865 1868
    }
1866 1869

	
1867 1870
    /// \brief Returns the value of the node.
1868 1871
    ///
1869 1872
    /// Returns the the value of the node.
1870 1873
    Value nodeValue(const Node& n) const {
1871 1874
      return (*_node_potential)[n];
1872 1875
    }
1873 1876

	
1874 1877
    /// \brief Returns the number of the blossoms in the basis.
1875 1878
    ///
1876 1879
    /// Returns the number of the blossoms in the basis.
1877 1880
    /// \see BlossomIt
1878 1881
    int blossomNum() const {
1879 1882
      return _blossom_potential.size();
1880 1883
    }
1881 1884

	
1882 1885

	
1883 1886
    /// \brief Returns the number of the nodes in the blossom.
1884 1887
    ///
1885 1888
    /// Returns the number of the nodes in the blossom.
1886 1889
    int blossomSize(int k) const {
1887 1890
      return _blossom_potential[k].end - _blossom_potential[k].begin;
1888 1891
    }
1889 1892

	
1890 1893
    /// \brief Returns the value of the blossom.
1891 1894
    ///
1892 1895
    /// Returns the the value of the blossom.
1893 1896
    /// \see BlossomIt
1894 1897
    Value blossomValue(int k) const {
1895 1898
      return _blossom_potential[k].value;
1896 1899
    }
1897 1900

	
1898 1901
    /// \brief Lemon iterator for get the items of the blossom.
1899 1902
    ///
1900 1903
    /// Lemon iterator for get the nodes of the blossom. This class
1901 1904
    /// provides a common style lemon iterator which gives back a
1902 1905
    /// subset of the nodes.
1903 1906
    class BlossomIt {
1904 1907
    public:
1905 1908

	
1906 1909
      /// \brief Constructor.
1907 1910
      ///
1908 1911
      /// Constructor for get the nodes of the variable.
1909 1912
      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
1910 1913
        : _algorithm(&algorithm)
1911 1914
      {
1912 1915
        _index = _algorithm->_blossom_potential[variable].begin;
1913 1916
        _last = _algorithm->_blossom_potential[variable].end;
1914 1917
      }
1915 1918

	
1916
      /// \brief Invalid constructor.
1917
      ///
1918
      /// Invalid constructor.
1919
      BlossomIt(Invalid) : _index(-1) {}
1920

	
1921 1919
      /// \brief Conversion to node.
1922 1920
      ///
1923 1921
      /// Conversion to node.
1924 1922
      operator Node() const {
1925
        return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
1923
        return _algorithm->_blossom_node_list[_index];
1926 1924
      }
1927 1925

	
1928 1926
      /// \brief Increment operator.
1929 1927
      ///
1930 1928
      /// Increment operator.
1931 1929
      BlossomIt& operator++() {
1932 1930
        ++_index;
1933
        if (_index == _last) {
1934
          _index = -1;
1935
        }
1936 1931
        return *this;
1937 1932
      }
1938 1933

	
1939
      bool operator==(const BlossomIt& it) const {
1940
        return _index == it._index;
1941
      }
1942
      bool operator!=(const BlossomIt& it) const {
1943
        return _index != it._index;
1944
      }
1934
      /// \brief Validity checking
1935
      ///
1936
      /// Checks whether the iterator is invalid.
1937
      bool operator==(Invalid) const { return _index == _last; }
1938

	
1939
      /// \brief Validity checking
1940
      ///
1941
      /// Checks whether the iterator is valid.
1942
      bool operator!=(Invalid) const { return _index != _last; }
1945 1943

	
1946 1944
    private:
1947 1945
      const MaxWeightedMatching* _algorithm;
1948 1946
      int _last;
1949 1947
      int _index;
1950 1948
    };
1951 1949

	
1952 1950
    /// @}
1953 1951

	
1954 1952
  };
1955 1953

	
1956 1954
  /// \ingroup matching
1957 1955
  ///
1958 1956
  /// \brief Weighted perfect matching in general graphs
1959 1957
  ///
1960 1958
  /// This class provides an efficient implementation of Edmond's
1961
  /// maximum weighted perfecr matching algorithm. The implementation
1959
  /// maximum weighted perfect matching algorithm. The implementation
1962 1960
  /// is based on extensive use of priority queues and provides
1963 1961
  /// \f$O(nm\log(n))\f$ time complexity.
1964 1962
  ///
1965 1963
  /// The maximum weighted matching problem is to find undirected
1966
  /// arcs in the digraph with maximum overall weight and no two of
1967
  /// them shares their endpoints and covers all nodes. The problem
1968
  /// can be formulated with the next linear program:
1964
  /// edges in the graph with maximum overall weight and no two of
1965
  /// them shares their ends and covers all nodes. The problem can be
1966
  /// formulated with the following linear program.
1969 1967
  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1970
  ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
1968
  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
1969
      \quad \forall B\in\mathcal{O}\f] */
1971 1970
  /// \f[x_e \ge 0\quad \forall e\in E\f]
1972 1971
  /// \f[\max \sum_{e\in E}x_ew_e\f]
1973
  /// where \f$\delta(X)\f$ is the set of arcs incident to a node in
1974
  /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in
1975
  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
1976
  /// the nodes.
1972
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1973
  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
1974
  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
1975
  /// subsets of the nodes.
1977 1976
  ///
1978 1977
  /// The algorithm calculates an optimal matching and a proof of the
1979 1978
  /// optimality. The solution of the dual problem can be used to check
1980
  /// the result of the algorithm. The dual linear problem is the next:
1981
  /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
1979
  /// the result of the algorithm. The dual linear problem is the
1980
  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
1981
      w_{uv} \quad \forall uv\in E\f] */
1982 1982
  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1983
  /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
1983
  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1984
      \frac{\vert B \vert - 1}{2}z_B\f] */
1984 1985
  ///
1985 1986
  /// The algorithm can be executed with \c run() or the \c init() and
1986 1987
  /// then the \c start() member functions. After it the matching can
1987 1988
  /// be asked with \c matching() or mate() functions. The dual
1988 1989
  /// solution can be get with \c nodeValue(), \c blossomNum() and \c
1989 1990
  /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
1990 1991
  /// "BlossomIt" nested class which is able to iterate on the nodes
1991 1992
  /// of a blossom. If the value type is integral then the dual
1992 1993
  /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
1993 1994
  template <typename _Graph,
1994 1995
            typename _WeightMap = typename _Graph::template EdgeMap<int> >
1995 1996
  class MaxWeightedPerfectMatching {
1996 1997
  public:
1997 1998

	
1998 1999
    typedef _Graph Graph;
1999 2000
    typedef _WeightMap WeightMap;
2000 2001
    typedef typename WeightMap::Value Value;
2001 2002

	
2002 2003
    /// \brief Scaling factor for dual solution
2003 2004
    ///
2004 2005
    /// Scaling factor for dual solution, it is equal to 4 or 1
2005 2006
    /// according to the value type.
2006 2007
    static const int dualScale =
2007 2008
      std::numeric_limits<Value>::is_integer ? 4 : 1;
2008 2009

	
2009 2010
    typedef typename Graph::template NodeMap<typename Graph::Arc>
2010 2011
    MatchingMap;
2011 2012

	
2012 2013
  private:
2013 2014

	
2014 2015
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
2015 2016

	
2016 2017
    typedef typename Graph::template NodeMap<Value> NodePotential;
2017 2018
    typedef std::vector<Node> BlossomNodeList;
2018 2019

	
2019 2020
    struct BlossomVariable {
2020 2021
      int begin, end;
2021 2022
      Value value;
2022 2023

	
2023 2024
      BlossomVariable(int _begin, int _end, Value _value)
2024 2025
        : begin(_begin), end(_end), value(_value) {}
2025 2026

	
2026 2027
    };
2027 2028

	
2028 2029
    typedef std::vector<BlossomVariable> BlossomPotential;
2029 2030

	
2030 2031
    const Graph& _graph;
2031 2032
    const WeightMap& _weight;
2032 2033

	
2033 2034
    MatchingMap* _matching;
2034 2035

	
2035 2036
    NodePotential* _node_potential;
2036 2037

	
2037 2038
    BlossomPotential _blossom_potential;
2038 2039
    BlossomNodeList _blossom_node_list;
2039 2040

	
2040 2041
    int _node_num;
2041 2042
    int _blossom_num;
2042 2043

	
2043
    typedef typename Graph::template NodeMap<int> NodeIntMap;
2044
    typedef typename Graph::template ArcMap<int> ArcIntMap;
2045
    typedef typename Graph::template EdgeMap<int> EdgeIntMap;
2046 2044
    typedef RangeMap<int> IntIntMap;
2047 2045

	
2048 2046
    enum Status {
2049 2047
      EVEN = -1, MATCHED = 0, ODD = 1
2050 2048
    };
2051 2049

	
2052
    typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
2050
    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
2053 2051
    struct BlossomData {
2054 2052
      int tree;
2055 2053
      Status status;
2056 2054
      Arc pred, next;
2057 2055
      Value pot, offset;
2058 2056
    };
2059 2057

	
2060
    NodeIntMap *_blossom_index;
2058
    IntNodeMap *_blossom_index;
2061 2059
    BlossomSet *_blossom_set;
2062 2060
    RangeMap<BlossomData>* _blossom_data;
2063 2061

	
2064
    NodeIntMap *_node_index;
2065
    ArcIntMap *_node_heap_index;
2062
    IntNodeMap *_node_index;
2063
    IntArcMap *_node_heap_index;
2066 2064

	
2067 2065
    struct NodeData {
2068 2066

	
2069
      NodeData(ArcIntMap& node_heap_index)
2067
      NodeData(IntArcMap& node_heap_index)
2070 2068
        : heap(node_heap_index) {}
2071 2069

	
2072 2070
      int blossom;
2073 2071
      Value pot;
2074
      BinHeap<Value, ArcIntMap> heap;
2072
      BinHeap<Value, IntArcMap> heap;
2075 2073
      std::map<int, Arc> heap_index;
2076 2074

	
2077 2075
      int tree;
2078 2076
    };
2079 2077

	
2080 2078
    RangeMap<NodeData>* _node_data;
2081 2079

	
2082 2080
    typedef ExtendFindEnum<IntIntMap> TreeSet;
2083 2081

	
2084 2082
    IntIntMap *_tree_set_index;
2085 2083
    TreeSet *_tree_set;
2086 2084

	
2087 2085
    IntIntMap *_delta2_index;
2088 2086
    BinHeap<Value, IntIntMap> *_delta2;
2089 2087

	
2090
    EdgeIntMap *_delta3_index;
2091
    BinHeap<Value, EdgeIntMap> *_delta3;
2088
    IntEdgeMap *_delta3_index;
2089
    BinHeap<Value, IntEdgeMap> *_delta3;
2092 2090

	
2093 2091
    IntIntMap *_delta4_index;
2094 2092
    BinHeap<Value, IntIntMap> *_delta4;
2095 2093

	
2096 2094
    Value _delta_sum;
2097 2095

	
2098 2096
    void createStructures() {
2099 2097
      _node_num = countNodes(_graph);
2100 2098
      _blossom_num = _node_num * 3 / 2;
2101 2099

	
2102 2100
      if (!_matching) {
2103 2101
        _matching = new MatchingMap(_graph);
2104 2102
      }
2105 2103
      if (!_node_potential) {
2106 2104
        _node_potential = new NodePotential(_graph);
2107 2105
      }
2108 2106
      if (!_blossom_set) {
2109
        _blossom_index = new NodeIntMap(_graph);
2107
        _blossom_index = new IntNodeMap(_graph);
2110 2108
        _blossom_set = new BlossomSet(*_blossom_index);
2111 2109
        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2112 2110
      }
2113 2111

	
2114 2112
      if (!_node_index) {
2115
        _node_index = new NodeIntMap(_graph);
2116
        _node_heap_index = new ArcIntMap(_graph);
2113
        _node_index = new IntNodeMap(_graph);
2114
        _node_heap_index = new IntArcMap(_graph);
2117 2115
        _node_data = new RangeMap<NodeData>(_node_num,
2118
                                              NodeData(*_node_heap_index));
2116
                                            NodeData(*_node_heap_index));
2119 2117
      }
2120 2118

	
2121 2119
      if (!_tree_set) {
2122 2120
        _tree_set_index = new IntIntMap(_blossom_num);
2123 2121
        _tree_set = new TreeSet(*_tree_set_index);
2124 2122
      }
2125 2123
      if (!_delta2) {
2126 2124
        _delta2_index = new IntIntMap(_blossom_num);
2127 2125
        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2128 2126
      }
2129 2127
      if (!_delta3) {
2130
        _delta3_index = new EdgeIntMap(_graph);
2131
        _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
2128
        _delta3_index = new IntEdgeMap(_graph);
2129
        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
2132 2130
      }
2133 2131
      if (!_delta4) {
2134 2132
        _delta4_index = new IntIntMap(_blossom_num);
2135 2133
        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2136 2134
      }
2137 2135
    }
2138 2136

	
2139 2137
    void destroyStructures() {
2140 2138
      _node_num = countNodes(_graph);
2141 2139
      _blossom_num = _node_num * 3 / 2;
2142 2140

	
2143 2141
      if (_matching) {
2144 2142
        delete _matching;
2145 2143
      }
2146 2144
      if (_node_potential) {
2147 2145
        delete _node_potential;
2148 2146
      }
2149 2147
      if (_blossom_set) {
2150 2148
        delete _blossom_index;
2151 2149
        delete _blossom_set;
2152 2150
        delete _blossom_data;
2153 2151
      }
2154 2152

	
2155 2153
      if (_node_index) {
2156 2154
        delete _node_index;
2157 2155
        delete _node_heap_index;
2158 2156
        delete _node_data;
2159 2157
      }
2160 2158

	
2161 2159
      if (_tree_set) {
2162 2160
        delete _tree_set_index;
2163 2161
        delete _tree_set;
2164 2162
      }
2165 2163
      if (_delta2) {
2166 2164
        delete _delta2_index;
2167 2165
        delete _delta2;
2168 2166
      }
2169 2167
      if (_delta3) {
2170 2168
        delete _delta3_index;
2171 2169
        delete _delta3;
2172 2170
      }
2173 2171
      if (_delta4) {
2174 2172
        delete _delta4_index;
2175 2173
        delete _delta4;
2176 2174
      }
2177 2175
    }
2178 2176

	
2179 2177
    void matchedToEven(int blossom, int tree) {
... ...
@@ -2416,131 +2414,131 @@
2416 2414
                               (*_blossom_data)[vb].offset);
2417 2415
                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2418 2416
                           (*_blossom_data)[vb].offset) {
2419 2417
                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2420 2418
                                   (*_blossom_data)[vb].offset);
2421 2419
                }
2422 2420
              }
2423 2421
            }
2424 2422
          }
2425 2423
        }
2426 2424
      }
2427 2425
      (*_blossom_data)[blossom].offset = 0;
2428 2426
    }
2429 2427

	
2430 2428
    void alternatePath(int even, int tree) {
2431 2429
      int odd;
2432 2430

	
2433 2431
      evenToMatched(even, tree);
2434 2432
      (*_blossom_data)[even].status = MATCHED;
2435 2433

	
2436 2434
      while ((*_blossom_data)[even].pred != INVALID) {
2437 2435
        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
2438 2436
        (*_blossom_data)[odd].status = MATCHED;
2439 2437
        oddToMatched(odd);
2440 2438
        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2441 2439

	
2442 2440
        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
2443 2441
        (*_blossom_data)[even].status = MATCHED;
2444 2442
        evenToMatched(even, tree);
2445 2443
        (*_blossom_data)[even].next =
2446 2444
          _graph.oppositeArc((*_blossom_data)[odd].pred);
2447 2445
      }
2448 2446

	
2449 2447
    }
2450 2448

	
2451 2449
    void destroyTree(int tree) {
2452 2450
      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2453 2451
        if ((*_blossom_data)[b].status == EVEN) {
2454 2452
          (*_blossom_data)[b].status = MATCHED;
2455 2453
          evenToMatched(b, tree);
2456 2454
        } else if ((*_blossom_data)[b].status == ODD) {
2457 2455
          (*_blossom_data)[b].status = MATCHED;
2458 2456
          oddToMatched(b);
2459 2457
        }
2460 2458
      }
2461 2459
      _tree_set->eraseClass(tree);
2462 2460
    }
2463 2461

	
2464
    void augmentOnArc(const Edge& arc) {
2465

	
2466
      int left = _blossom_set->find(_graph.u(arc));
2467
      int right = _blossom_set->find(_graph.v(arc));
2462
    void augmentOnEdge(const Edge& edge) {
2463

	
2464
      int left = _blossom_set->find(_graph.u(edge));
2465
      int right = _blossom_set->find(_graph.v(edge));
2468 2466

	
2469 2467
      int left_tree = _tree_set->find(left);
2470 2468
      alternatePath(left, left_tree);
2471 2469
      destroyTree(left_tree);
2472 2470

	
2473 2471
      int right_tree = _tree_set->find(right);
2474 2472
      alternatePath(right, right_tree);
2475 2473
      destroyTree(right_tree);
2476 2474

	
2477
      (*_blossom_data)[left].next = _graph.direct(arc, true);
2478
      (*_blossom_data)[right].next = _graph.direct(arc, false);
2475
      (*_blossom_data)[left].next = _graph.direct(edge, true);
2476
      (*_blossom_data)[right].next = _graph.direct(edge, false);
2479 2477
    }
2480 2478

	
2481 2479
    void extendOnArc(const Arc& arc) {
2482 2480
      int base = _blossom_set->find(_graph.target(arc));
2483 2481
      int tree = _tree_set->find(base);
2484 2482

	
2485 2483
      int odd = _blossom_set->find(_graph.source(arc));
2486 2484
      _tree_set->insert(odd, tree);
2487 2485
      (*_blossom_data)[odd].status = ODD;
2488 2486
      matchedToOdd(odd);
2489 2487
      (*_blossom_data)[odd].pred = arc;
2490 2488

	
2491 2489
      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
2492 2490
      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2493 2491
      _tree_set->insert(even, tree);
2494 2492
      (*_blossom_data)[even].status = EVEN;
2495 2493
      matchedToEven(even, tree);
2496 2494
    }
2497 2495

	
2498
    void shrinkOnArc(const Edge& edge, int tree) {
2496
    void shrinkOnEdge(const Edge& edge, int tree) {
2499 2497
      int nca = -1;
2500 2498
      std::vector<int> left_path, right_path;
2501 2499

	
2502 2500
      {
2503 2501
        std::set<int> left_set, right_set;
2504 2502
        int left = _blossom_set->find(_graph.u(edge));
2505 2503
        left_path.push_back(left);
2506 2504
        left_set.insert(left);
2507 2505

	
2508 2506
        int right = _blossom_set->find(_graph.v(edge));
2509 2507
        right_path.push_back(right);
2510 2508
        right_set.insert(right);
2511 2509

	
2512 2510
        while (true) {
2513 2511

	
2514 2512
          if ((*_blossom_data)[left].pred == INVALID) break;
2515 2513

	
2516 2514
          left =
2517 2515
            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2518 2516
          left_path.push_back(left);
2519 2517
          left =
2520 2518
            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2521 2519
          left_path.push_back(left);
2522 2520

	
2523 2521
          left_set.insert(left);
2524 2522

	
2525 2523
          if (right_set.find(left) != right_set.end()) {
2526 2524
            nca = left;
2527 2525
            break;
2528 2526
          }
2529 2527

	
2530 2528
          if ((*_blossom_data)[right].pred == INVALID) break;
2531 2529

	
2532 2530
          right =
2533 2531
            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2534 2532
          right_path.push_back(right);
2535 2533
          right =
2536 2534
            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2537 2535
          right_path.push_back(right);
2538 2536

	
2539 2537
          right_set.insert(right);
2540 2538

	
2541 2539
          if (left_set.find(right) != left_set.end()) {
2542 2540
            nca = right;
2543 2541
            break;
2544 2542
          }
2545 2543

	
2546 2544
        }
... ...
@@ -2786,327 +2784,322 @@
2786 2784
        for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2787 2785
             n != INVALID; ++n) {
2788 2786
          (*_node_data)[(*_node_index)[n]].pot -= offset;
2789 2787
        }
2790 2788

	
2791 2789
        Arc matching = (*_blossom_data)[blossoms[i]].next;
2792 2790
        Node base = _graph.source(matching);
2793 2791
        extractBlossom(blossoms[i], base, matching);
2794 2792
      }
2795 2793
    }
2796 2794

	
2797 2795
  public:
2798 2796

	
2799 2797
    /// \brief Constructor
2800 2798
    ///
2801 2799
    /// Constructor.
2802 2800
    MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2803 2801
      : _graph(graph), _weight(weight), _matching(0),
2804 2802
        _node_potential(0), _blossom_potential(), _blossom_node_list(),
2805 2803
        _node_num(0), _blossom_num(0),
2806 2804

	
2807 2805
        _blossom_index(0), _blossom_set(0), _blossom_data(0),
2808 2806
        _node_index(0), _node_heap_index(0), _node_data(0),
2809 2807
        _tree_set_index(0), _tree_set(0),
2810 2808

	
2811 2809
        _delta2_index(0), _delta2(0),
2812 2810
        _delta3_index(0), _delta3(0),
2813 2811
        _delta4_index(0), _delta4(0),
2814 2812

	
2815 2813
        _delta_sum() {}
2816 2814

	
2817 2815
    ~MaxWeightedPerfectMatching() {
2818 2816
      destroyStructures();
2819 2817
    }
2820 2818

	
2821 2819
    /// \name Execution control
2822 2820
    /// The simplest way to execute the algorithm is to use the member
2823 2821
    /// \c run() member function.
2824 2822

	
2825 2823
    ///@{
2826 2824

	
2827 2825
    /// \brief Initialize the algorithm
2828 2826
    ///
2829 2827
    /// Initialize the algorithm
2830 2828
    void init() {
2831 2829
      createStructures();
2832 2830

	
2833 2831
      for (ArcIt e(_graph); e != INVALID; ++e) {
2834
        _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
2832
        _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
2835 2833
      }
2836 2834
      for (EdgeIt e(_graph); e != INVALID; ++e) {
2837 2835
        _delta3_index->set(e, _delta3->PRE_HEAP);
2838 2836
      }
2839 2837
      for (int i = 0; i < _blossom_num; ++i) {
2840 2838
        _delta2_index->set(i, _delta2->PRE_HEAP);
2841 2839
        _delta4_index->set(i, _delta4->PRE_HEAP);
2842 2840
      }
2843 2841

	
2844 2842
      int index = 0;
2845 2843
      for (NodeIt n(_graph); n != INVALID; ++n) {
2846 2844
        Value max = - std::numeric_limits<Value>::max();
2847 2845
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2848 2846
          if (_graph.target(e) == n) continue;
2849 2847
          if ((dualScale * _weight[e]) / 2 > max) {
2850 2848
            max = (dualScale * _weight[e]) / 2;
2851 2849
          }
2852 2850
        }
2853 2851
        _node_index->set(n, index);
2854 2852
        (*_node_data)[index].pot = max;
2855 2853
        int blossom =
2856 2854
          _blossom_set->insert(n, std::numeric_limits<Value>::max());
2857 2855

	
2858 2856
        _tree_set->insert(blossom);
2859 2857

	
2860 2858
        (*_blossom_data)[blossom].status = EVEN;
2861 2859
        (*_blossom_data)[blossom].pred = INVALID;
2862 2860
        (*_blossom_data)[blossom].next = INVALID;
2863 2861
        (*_blossom_data)[blossom].pot = 0;
2864 2862
        (*_blossom_data)[blossom].offset = 0;
2865 2863
        ++index;
2866 2864
      }
2867 2865
      for (EdgeIt e(_graph); e != INVALID; ++e) {
2868 2866
        int si = (*_node_index)[_graph.u(e)];
2869 2867
        int ti = (*_node_index)[_graph.v(e)];
2870 2868
        if (_graph.u(e) != _graph.v(e)) {
2871 2869
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
2872 2870
                            dualScale * _weight[e]) / 2);
2873 2871
        }
2874 2872
      }
2875 2873
    }
2876 2874

	
2877 2875
    /// \brief Starts the algorithm
2878 2876
    ///
2879 2877
    /// Starts the algorithm
2880 2878
    bool start() {
2881 2879
      enum OpType {
2882 2880
        D2, D3, D4
2883 2881
      };
2884 2882

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

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

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

	
2896 2894
        _delta_sum = d2; OpType ot = D2;
2897 2895
        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
2898 2896
        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
2899 2897

	
2900 2898
        if (_delta_sum == std::numeric_limits<Value>::max()) {
2901 2899
          return false;
2902 2900
        }
2903 2901

	
2904 2902
        switch (ot) {
2905 2903
        case D2:
2906 2904
          {
2907 2905
            int blossom = _delta2->top();
2908 2906
            Node n = _blossom_set->classTop(blossom);
2909 2907
            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
2910 2908
            extendOnArc(e);
2911 2909
          }
2912 2910
          break;
2913 2911
        case D3:
2914 2912
          {
2915 2913
            Edge e = _delta3->top();
2916 2914

	
2917 2915
            int left_blossom = _blossom_set->find(_graph.u(e));
2918 2916
            int right_blossom = _blossom_set->find(_graph.v(e));
2919 2917

	
2920 2918
            if (left_blossom == right_blossom) {
2921 2919
              _delta3->pop();
2922 2920
            } else {
2923 2921
              int left_tree = _tree_set->find(left_blossom);
2924 2922
              int right_tree = _tree_set->find(right_blossom);
2925 2923

	
2926 2924
              if (left_tree == right_tree) {
2927
                shrinkOnArc(e, left_tree);
2925
                shrinkOnEdge(e, left_tree);
2928 2926
              } else {
2929
                augmentOnArc(e);
2927
                augmentOnEdge(e);
2930 2928
                unmatched -= 2;
2931 2929
              }
2932 2930
            }
2933 2931
          } break;
2934 2932
        case D4:
2935 2933
          splitBlossom(_delta4->top());
2936 2934
          break;
2937 2935
        }
2938 2936
      }
2939 2937
      extractMatching();
2940 2938
      return true;
2941 2939
    }
2942 2940

	
2943 2941
    /// \brief Runs %MaxWeightedPerfectMatching algorithm.
2944 2942
    ///
2945 2943
    /// This method runs the %MaxWeightedPerfectMatching algorithm.
2946 2944
    ///
2947 2945
    /// \note mwm.run() is just a shortcut of the following code.
2948 2946
    /// \code
2949 2947
    ///   mwm.init();
2950 2948
    ///   mwm.start();
2951 2949
    /// \endcode
2952 2950
    bool run() {
2953 2951
      init();
2954 2952
      return start();
2955 2953
    }
2956 2954

	
2957 2955
    /// @}
2958 2956

	
2959 2957
    /// \name Primal solution
2960 2958
    /// Functions for get the primal solution, ie. the matching.
2961 2959

	
2962 2960
    /// @{
2963 2961

	
2964 2962
    /// \brief Returns the matching value.
2965 2963
    ///
2966 2964
    /// Returns the matching value.
2967 2965
    Value matchingValue() const {
2968 2966
      Value sum = 0;
2969 2967
      for (NodeIt n(_graph); n != INVALID; ++n) {
2970 2968
        if ((*_matching)[n] != INVALID) {
2971 2969
          sum += _weight[(*_matching)[n]];
2972 2970
        }
2973 2971
      }
2974 2972
      return sum /= 2;
2975 2973
    }
2976 2974

	
2977
    /// \brief Returns true when the arc is in the matching.
2975
    /// \brief Returns true when the edge is in the matching.
2978 2976
    ///
2979
    /// Returns true when the arc is in the matching.
2980
    bool matching(const Edge& arc) const {
2981
      return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
2977
    /// Returns true when the edge is in the matching.
2978
    bool matching(const Edge& edge) const {
2979
      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
2982 2980
    }
2983 2981

	
2984
    /// \brief Returns the incident matching arc.
2982
    /// \brief Returns the incident matching edge.
2985 2983
    ///
2986
    /// Returns the incident matching arc from given node.
2984
    /// Returns the incident matching arc from given edge.
2987 2985
    Arc matching(const Node& node) const {
2988 2986
      return (*_matching)[node];
2989 2987
    }
2990 2988

	
2991 2989
    /// \brief Returns the mate of the node.
2992 2990
    ///
2993 2991
    /// Returns the adjancent node in a mathcing arc.
2994 2992
    Node mate(const Node& node) const {
2995 2993
      return _graph.target((*_matching)[node]);
2996 2994
    }
2997 2995

	
2998 2996
    /// @}
2999 2997

	
3000 2998
    /// \name Dual solution
3001 2999
    /// Functions for get the dual solution.
3002 3000

	
3003 3001
    /// @{
3004 3002

	
3005 3003
    /// \brief Returns the value of the dual solution.
3006 3004
    ///
3007 3005
    /// Returns the value of the dual solution. It should be equal to
3008 3006
    /// the primal value scaled by \ref dualScale "dual scale".
3009 3007
    Value dualValue() const {
3010 3008
      Value sum = 0;
3011 3009
      for (NodeIt n(_graph); n != INVALID; ++n) {
3012 3010
        sum += nodeValue(n);
3013 3011
      }
3014 3012
      for (int i = 0; i < blossomNum(); ++i) {
3015 3013
        sum += blossomValue(i) * (blossomSize(i) / 2);
3016 3014
      }
3017 3015
      return sum;
3018 3016
    }
3019 3017

	
3020 3018
    /// \brief Returns the value of the node.
3021 3019
    ///
3022 3020
    /// Returns the the value of the node.
3023 3021
    Value nodeValue(const Node& n) const {
3024 3022
      return (*_node_potential)[n];
3025 3023
    }
3026 3024

	
3027 3025
    /// \brief Returns the number of the blossoms in the basis.
3028 3026
    ///
3029 3027
    /// Returns the number of the blossoms in the basis.
3030 3028
    /// \see BlossomIt
3031 3029
    int blossomNum() const {
3032 3030
      return _blossom_potential.size();
3033 3031
    }
3034 3032

	
3035 3033

	
3036 3034
    /// \brief Returns the number of the nodes in the blossom.
3037 3035
    ///
3038 3036
    /// Returns the number of the nodes in the blossom.
3039 3037
    int blossomSize(int k) const {
3040 3038
      return _blossom_potential[k].end - _blossom_potential[k].begin;
3041 3039
    }
3042 3040

	
3043 3041
    /// \brief Returns the value of the blossom.
3044 3042
    ///
3045 3043
    /// Returns the the value of the blossom.
3046 3044
    /// \see BlossomIt
3047 3045
    Value blossomValue(int k) const {
3048 3046
      return _blossom_potential[k].value;
3049 3047
    }
3050 3048

	
3051 3049
    /// \brief Lemon iterator for get the items of the blossom.
3052 3050
    ///
3053 3051
    /// Lemon iterator for get the nodes of the blossom. This class
3054 3052
    /// provides a common style lemon iterator which gives back a
3055 3053
    /// subset of the nodes.
3056 3054
    class BlossomIt {
3057 3055
    public:
3058 3056

	
3059 3057
      /// \brief Constructor.
3060 3058
      ///
3061 3059
      /// Constructor for get the nodes of the variable.
3062 3060
      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
3063 3061
        : _algorithm(&algorithm)
3064 3062
      {
3065 3063
        _index = _algorithm->_blossom_potential[variable].begin;
3066 3064
        _last = _algorithm->_blossom_potential[variable].end;
3067 3065
      }
3068 3066

	
3069
      /// \brief Invalid constructor.
3070
      ///
3071
      /// Invalid constructor.
3072
      BlossomIt(Invalid) : _index(-1) {}
3073

	
3074 3067
      /// \brief Conversion to node.
3075 3068
      ///
3076 3069
      /// Conversion to node.
3077 3070
      operator Node() const {
3078
        return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
3071
        return _algorithm->_blossom_node_list[_index];
3079 3072
      }
3080 3073

	
3081 3074
      /// \brief Increment operator.
3082 3075
      ///
3083 3076
      /// Increment operator.
3084 3077
      BlossomIt& operator++() {
3085 3078
        ++_index;
3086
        if (_index == _last) {
3087
          _index = -1;
3088
        }
3089 3079
        return *this;
3090 3080
      }
3091 3081

	
3092
      bool operator==(const BlossomIt& it) const {
3093
        return _index == it._index;
3094
      }
3095
      bool operator!=(const BlossomIt& it) const {
3096
        return _index != it._index;
3097
      }
3082
      /// \brief Validity checking
3083
      ///
3084
      /// Checks whether the iterator is invalid.
3085
      bool operator==(Invalid) const { return _index == _last; }
3086

	
3087
      /// \brief Validity checking
3088
      ///
3089
      /// Checks whether the iterator is valid.
3090
      bool operator!=(Invalid) const { return _index != _last; }
3098 3091

	
3099 3092
    private:
3100 3093
      const MaxWeightedPerfectMatching* _algorithm;
3101 3094
      int _last;
3102 3095
      int _index;
3103 3096
    };
3104 3097

	
3105 3098
    /// @}
3106 3099

	
3107 3100
  };
3108 3101

	
3109 3102

	
3110 3103
} //END OF NAMESPACE LEMON
3111 3104

	
3112 3105
#endif //LEMON_MAX_MATCHING_H
Ignore white space 96 line context
1 1
INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR})
2 2

	
3 3
LINK_DIRECTORIES(${CMAKE_BINARY_DIR}/lemon)
4 4

	
5 5
SET(TESTS
6 6
  bfs_test
7 7
  counter_test
8 8
  dfs_test
9 9
  digraph_test
10 10
  dijkstra_test
11 11
  dim_test
12 12
  error_test
13 13
  graph_copy_test
14 14
  graph_test
15 15
  graph_utils_test
16 16
  heap_test
17 17
  kruskal_test
18 18
  maps_test
19 19
  max_matching_test
20
  max_weighted_matching_test
21 20
  random_test
22 21
  path_test
23 22
  time_measure_test
24 23
  unionfind_test)
25 24

	
26 25
FOREACH(TEST_NAME ${TESTS})
27 26
  ADD_EXECUTABLE(${TEST_NAME} ${TEST_NAME}.cc)
28 27
  TARGET_LINK_LIBRARIES(${TEST_NAME} lemon)
29 28
  ADD_TEST(${TEST_NAME} ${TEST_NAME})
30 29
ENDFOREACH(TEST_NAME)
Ignore white space 6 line context
1 1
EXTRA_DIST += \
2 2
	test/CMakeLists.txt
3 3

	
4 4
noinst_HEADERS += \
5 5
	test/graph_test.h \
6 6
        test/test_tools.h
7 7

	
8 8
check_PROGRAMS += \
9 9
	test/bfs_test \
10 10
        test/counter_test \
11 11
	test/dfs_test \
12 12
	test/digraph_test \
13 13
	test/dijkstra_test \
14 14
        test/dim_test \
15 15
	test/error_test \
16 16
	test/graph_copy_test \
17 17
	test/graph_test \
18 18
	test/graph_utils_test \
19 19
	test/heap_test \
20 20
	test/kruskal_test \
21 21
        test/maps_test \
22 22
	test/max_matching_test \
23
	test/max_weighted_matching_test \
24 23
        test/random_test \
25 24
        test/path_test \
26 25
        test/test_tools_fail \
27 26
        test/test_tools_pass \
28 27
        test/time_measure_test \
29 28
	test/unionfind_test
30 29

	
31 30
TESTS += $(check_PROGRAMS)
32 31
XFAIL_TESTS += test/test_tools_fail$(EXEEXT)
33 32

	
34 33
test_bfs_test_SOURCES = test/bfs_test.cc
35 34
test_counter_test_SOURCES = test/counter_test.cc
36 35
test_dfs_test_SOURCES = test/dfs_test.cc
37 36
test_digraph_test_SOURCES = test/digraph_test.cc
38 37
test_dijkstra_test_SOURCES = test/dijkstra_test.cc
39 38
test_dim_test_SOURCES = test/dim_test.cc
40 39
test_error_test_SOURCES = test/error_test.cc
41 40
test_graph_copy_test_SOURCES = test/graph_copy_test.cc
42 41
test_graph_test_SOURCES = test/graph_test.cc
43 42
test_graph_utils_test_SOURCES = test/graph_utils_test.cc
44 43
test_heap_test_SOURCES = test/heap_test.cc
45 44
test_kruskal_test_SOURCES = test/kruskal_test.cc
46 45
test_maps_test_SOURCES = test/maps_test.cc
47 46
test_max_matching_test_SOURCES = test/max_matching_test.cc
48
test_max_weighted_matching_test_SOURCES = test/max_weighted_matching_test.cc
49 47
test_path_test_SOURCES = test/path_test.cc
50 48
test_random_test_SOURCES = test/random_test.cc
51 49
test_test_tools_fail_SOURCES = test/test_tools_fail.cc
52 50
test_test_tools_pass_SOURCES = test/test_tools_pass.cc
53 51
test_time_measure_test_SOURCES = test/time_measure_test.cc
54 52
test_unionfind_test_SOURCES = test/unionfind_test.cc
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-2008
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
#include <sstream>
20 21
#include <vector>
21 22
#include <queue>
22 23
#include <lemon/math.h>
23 24
#include <cstdlib>
24 25

	
26
#include <lemon/max_matching.h>
27
#include <lemon/smart_graph.h>
28
#include <lemon/lgf_reader.h>
29

	
25 30
#include "test_tools.h"
26
#include <lemon/list_graph.h>
27
#include <lemon/max_matching.h>
28 31

	
29 32
using namespace std;
30 33
using namespace lemon;
31 34

	
35
GRAPH_TYPEDEFS(SmartGraph);
36

	
37

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

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

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

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

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

	
120
  int barrier_num = 0;
121

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

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

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

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

	
152
  std::set<int> comp_root;
153
  int odd_comp_num = 0;
154
  for (NodeIt n(graph); n != INVALID; ++n) {
155
    if (mm.decomposition(n) != MaxMatching<SmartGraph>::ODD) {
156
      int root = comp.find(n);
157
      if (comp_root.find(root) == comp_root.end()) {
158
        comp_root.insert(root);
159
        if (comp.size(n) % 2 == 1) {
160
          ++odd_comp_num;
161
        }
162
      }
163
    }
164
  }
165

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

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

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

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

	
197
  int pv = 0;
198
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
199
    if (mwm.matching(n) != INVALID) {
200
      check(mwm.nodeValue(n) >= 0, "Invalid node value");
201
      pv += weight[mwm.matching(n)];
202
      SmartGraph::Node o = graph.target(mwm.matching(n));
203
      check(mwm.mate(n) == o, "Invalid matching");
204
      check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
205
            "Invalid matching");
206
    } else {
207
      check(mwm.mate(n) == INVALID, "Invalid matching");
208
      check(mwm.nodeValue(n) == 0, "Invalid matching");
209
    }
210
  }
211

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

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

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

	
225
  return;
226
}
227

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

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

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

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

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

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

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

	
276
  return;
277
}
278

	
279

	
32 280
int main() {
33 281

	
34
  typedef ListGraph Graph;
282
  for (int i = 0; i < lgfn; ++i) {
283
    SmartGraph graph;
284
    SmartGraph::EdgeMap<int> weight(graph);
35 285

	
36
  GRAPH_TYPEDEFS(Graph);
286
    istringstream lgfs(lgf[i]);
287
    graphReader(graph, lgfs).
288
      edgeMap("weight", weight).run();
37 289

	
38
  Graph g;
39
  g.clear();
290
    MaxMatching<SmartGraph> mm(graph);
291
    mm.run();
292
    checkMatching(graph, mm);
40 293

	
41
  std::vector<Graph::Node> nodes;
42
  for (int i=0; i<13; ++i)
43
      nodes.push_back(g.addNode());
294
    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
295
    mwm.run();
296
    checkWeightedMatching(graph, weight, mwm);
44 297

	
45
  g.addEdge(nodes[0], nodes[0]);
46
  g.addEdge(nodes[6], nodes[10]);
47
  g.addEdge(nodes[5], nodes[10]);
48
  g.addEdge(nodes[4], nodes[10]);
49
  g.addEdge(nodes[3], nodes[11]);
50
  g.addEdge(nodes[1], nodes[6]);
51
  g.addEdge(nodes[4], nodes[7]);
52
  g.addEdge(nodes[1], nodes[8]);
53
  g.addEdge(nodes[0], nodes[8]);
54
  g.addEdge(nodes[3], nodes[12]);
55
  g.addEdge(nodes[6], nodes[9]);
56
  g.addEdge(nodes[9], nodes[11]);
57
  g.addEdge(nodes[2], nodes[10]);
58
  g.addEdge(nodes[10], nodes[8]);
59
  g.addEdge(nodes[5], nodes[8]);
60
  g.addEdge(nodes[6], nodes[3]);
61
  g.addEdge(nodes[0], nodes[5]);
62
  g.addEdge(nodes[6], nodes[12]);
298
    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
299
    bool perfect = mwpm.run();
63 300

	
64
  MaxMatching<Graph> max_matching(g);
65
  max_matching.init();
66
  max_matching.startDense();
301
    check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
302
          "Perfect matching found");
67 303

	
68
  int s=0;
69
  Graph::NodeMap<Node> mate(g,INVALID);
70
  max_matching.mateMap(mate);
71
  for(NodeIt v(g); v!=INVALID; ++v) {
72
    if ( mate[v]!=INVALID ) ++s;
73
  }
74
  int size=int(s/2);  //size will be used as the size of a maxmatching
75

	
76
  for(NodeIt v(g); v!=INVALID; ++v) {
77
    max_matching.mate(v);
78
  }
79

	
80
  check ( size == max_matching.size(), "mate() returns a different size matching than max_matching.size()" );
81

	
82
  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos0(g);
83
  max_matching.decomposition(pos0);
84

	
85
  max_matching.init();
86
  max_matching.startSparse();
87
  s=0;
88
  max_matching.mateMap(mate);
89
  for(NodeIt v(g); v!=INVALID; ++v) {
90
    if ( mate[v]!=INVALID ) ++s;
91
  }
92
  check ( int(s/2) == size, "The size does not equal!" );
93

	
94
  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos1(g);
95
  max_matching.decomposition(pos1);
96

	
97
  max_matching.run();
98
  s=0;
99
  max_matching.mateMap(mate);
100
  for(NodeIt v(g); v!=INVALID; ++v) {
101
    if ( mate[v]!=INVALID ) ++s;
102
  }
103
  check ( int(s/2) == size, "The size does not equal!" );
104

	
105
  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos2(g);
106
  max_matching.decomposition(pos2);
107

	
108
  max_matching.run();
109
  s=0;
110
  max_matching.mateMap(mate);
111
  for(NodeIt v(g); v!=INVALID; ++v) {
112
    if ( mate[v]!=INVALID ) ++s;
113
  }
114
  check ( int(s/2) == size, "The size does not equal!" );
115

	
116
  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos(g);
117
  max_matching.decomposition(pos);
118

	
119
  bool ismatching=true;
120
  for(NodeIt v(g); v!=INVALID; ++v) {
121
    if ( mate[v]!=INVALID ) {
122
      Node u=mate[v];
123
      if (mate[u]!=v) ismatching=false;
304
    if (perfect) {
305
      checkWeightedPerfectMatching(graph, weight, mwpm);
124 306
    }
125 307
  }
126
  check ( ismatching, "It is not a matching!" );
127

	
128
  bool coincide=true;
129
  for(NodeIt v(g); v!=INVALID; ++v) {
130
   if ( pos0[v] != pos1[v] || pos1[v]!=pos2[v] || pos2[v]!=pos[v] ) {
131
     coincide=false;
132
    }
133
  }
134
  check ( coincide, "The decompositions do not coincide! " );
135

	
136
  bool noarc=true;
137
  for(EdgeIt e(g); e!=INVALID; ++e) {
138
   if ( (pos[g.v(e)]==max_matching.C &&
139
         pos[g.u(e)]==max_matching.D) ||
140
         (pos[g.v(e)]==max_matching.D &&
141
          pos[g.u(e)]==max_matching.C) )
142
      noarc=false;
143
  }
144
  check ( noarc, "There are arcs between D and C!" );
145

	
146
  bool oddcomp=true;
147
  Graph::NodeMap<bool> todo(g,true);
148
  int num_comp=0;
149
  for(NodeIt v(g); v!=INVALID; ++v) {
150
   if ( pos[v]==max_matching.D && todo[v] ) {
151
      int comp_size=1;
152
      ++num_comp;
153
      std::queue<Node> Q;
154
      Q.push(v);
155
      todo.set(v,false);
156
      while (!Q.empty()) {
157
        Node w=Q.front();
158
        Q.pop();
159
        for(IncEdgeIt e(g,w); e!=INVALID; ++e) {
160
          Node u=g.runningNode(e);
161
          if ( pos[u]==max_matching.D && todo[u] ) {
162
            ++comp_size;
163
            Q.push(u);
164
            todo.set(u,false);
165
          }
166
        }
167
      }
168
      if ( !(comp_size % 2) ) oddcomp=false;
169
    }
170
  }
171
  check ( oddcomp, "A component of g[D] is not odd." );
172

	
173
  int barrier=0;
174
  for(NodeIt v(g); v!=INVALID; ++v) {
175
    if ( pos[v]==max_matching.A ) ++barrier;
176
  }
177
  int expected_size=int( countNodes(g)-num_comp+barrier)/2;
178
  check ( size==expected_size, "The size of the matching is wrong." );
179 308

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

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

	
26
#include "test_tools.h"
27
#include <lemon/smart_graph.h>
28
#include <lemon/max_matching.h>
29
#include <lemon/lgf_reader.h>
30

	
31
using namespace std;
32
using namespace lemon;
33

	
34
GRAPH_TYPEDEFS(SmartGraph);
35

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

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

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

	
111
void checkMatching(const SmartGraph& graph,
112
                   const SmartGraph::EdgeMap<int>& weight,
113
                   const MaxWeightedMatching<SmartGraph>& mwm) {
114
  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
115
    if (graph.u(e) == graph.v(e)) continue;
116
    int rw =
117
      mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
118

	
119
    for (int i = 0; i < mwm.blossomNum(); ++i) {
120
      bool s = false, t = false;
121
      for (MaxWeightedMatching<SmartGraph>::BlossomIt n(mwm, i);
122
           n != INVALID; ++n) {
123
        if (graph.u(e) == n) s = true;
124
        if (graph.v(e) == n) t = true;
125
      }
126
      if (s == true && t == true) {
127
        rw += mwm.blossomValue(i);
128
      }
129
    }
130
    rw -= weight[e] * mwm.dualScale;
131

	
132
    check(rw >= 0, "Negative reduced weight");
133
    check(rw == 0 || !mwm.matching(e),
134
          "Non-zero reduced weight on matching arc");
135
  }
136

	
137
  int pv = 0;
138
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
139
    if (mwm.matching(n) != INVALID) {
140
      check(mwm.nodeValue(n) >= 0, "Invalid node value");
141
      pv += weight[mwm.matching(n)];
142
      SmartGraph::Node o = graph.target(mwm.matching(n));
143
      check(mwm.mate(n) == o, "Invalid matching");
144
      check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
145
            "Invalid matching");
146
    } else {
147
      check(mwm.mate(n) == INVALID, "Invalid matching");
148
      check(mwm.nodeValue(n) == 0, "Invalid matching");
149
    }
150
  }
151

	
152
  int dv = 0;
153
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
154
    dv += mwm.nodeValue(n);
155
  }
156

	
157
  for (int i = 0; i < mwm.blossomNum(); ++i) {
158
    check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
159
    check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
160
    dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
161
  }
162

	
163
  check(pv * mwm.dualScale == dv * 2, "Wrong duality");
164

	
165
  return;
166
}
167

	
168
void checkPerfectMatching(const SmartGraph& graph,
169
                          const SmartGraph::EdgeMap<int>& weight,
170
                          const MaxWeightedPerfectMatching<SmartGraph>& mwpm) {
171
  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
172
    if (graph.u(e) == graph.v(e)) continue;
173
    int rw =
174
      mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
175

	
176
    for (int i = 0; i < mwpm.blossomNum(); ++i) {
177
      bool s = false, t = false;
178
      for (MaxWeightedPerfectMatching<SmartGraph>::BlossomIt n(mwpm, i);
179
           n != INVALID; ++n) {
180
        if (graph.u(e) == n) s = true;
181
        if (graph.v(e) == n) t = true;
182
      }
183
      if (s == true && t == true) {
184
        rw += mwpm.blossomValue(i);
185
      }
186
    }
187
    rw -= weight[e] * mwpm.dualScale;
188

	
189
    check(rw >= 0, "Negative reduced weight");
190
    check(rw == 0 || !mwpm.matching(e),
191
          "Non-zero reduced weight on matching arc");
192
  }
193

	
194
  int pv = 0;
195
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
196
    check(mwpm.matching(n) != INVALID, "Non perfect");
197
    pv += weight[mwpm.matching(n)];
198
    SmartGraph::Node o = graph.target(mwpm.matching(n));
199
    check(mwpm.mate(n) == o, "Invalid matching");
200
    check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
201
          "Invalid matching");
202
  }
203

	
204
  int dv = 0;
205
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
206
    dv += mwpm.nodeValue(n);
207
  }
208

	
209
  for (int i = 0; i < mwpm.blossomNum(); ++i) {
210
    check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
211
    check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
212
    dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
213
  }
214

	
215
  check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
216

	
217
  return;
218
}
219

	
220

	
221
int main() {
222

	
223
  for (int i = 0; i < lgfn; ++i) {
224
    SmartGraph graph;
225
    SmartGraph::EdgeMap<int> weight(graph);
226

	
227
    istringstream lgfs(lgf[i]);
228
    graphReader(graph, lgfs).
229
      edgeMap("weight", weight).run();
230

	
231
    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
232
    mwm.run();
233
    checkMatching(graph, weight, mwm);
234

	
235
    MaxMatching<SmartGraph> mm(graph);
236
    mm.run();
237

	
238
    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
239
    bool perfect = mwpm.run();
240

	
241
    check(perfect == (mm.size() * 2 == countNodes(graph)),
242
          "Perfect matching found");
243

	
244
    if (perfect) {
245
      checkPerfectMatching(graph, weight, mwpm);
246
    }
247
  }
248

	
249
  return 0;
250
}
0 comments (0 inline)