lemon/max_matching.h
author Balazs Dezso <deba@inf.elte.hu>
Sat, 18 Apr 2009 21:01:12 +0200
changeset 597 2ca0cdb5f366
parent 559 c5fd2d996909
child 590 b61354458b59
permissions -rw-r--r--
Fix in HaoOrlin (#264)
deba@326
     1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
deba@326
     2
 *
deba@326
     3
 * This file is a part of LEMON, a generic C++ optimization library.
deba@326
     4
 *
alpar@440
     5
 * Copyright (C) 2003-2009
deba@326
     6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
deba@326
     7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
deba@326
     8
 *
deba@326
     9
 * Permission to use, modify and distribute this software is granted
deba@326
    10
 * provided that this copyright notice appears in all copies. For
deba@326
    11
 * precise terms see the accompanying LICENSE file.
deba@326
    12
 *
deba@326
    13
 * This software is provided "AS IS" with no warranty of any kind,
deba@326
    14
 * express or implied, and with no claim as to its suitability for any
deba@326
    15
 * purpose.
deba@326
    16
 *
deba@326
    17
 */
deba@326
    18
deba@326
    19
#ifndef LEMON_MAX_MATCHING_H
deba@326
    20
#define LEMON_MAX_MATCHING_H
deba@326
    21
deba@326
    22
#include <vector>
deba@326
    23
#include <queue>
deba@326
    24
#include <set>
deba@326
    25
#include <limits>
deba@326
    26
deba@326
    27
#include <lemon/core.h>
deba@326
    28
#include <lemon/unionfind.h>
deba@326
    29
#include <lemon/bin_heap.h>
deba@326
    30
#include <lemon/maps.h>
deba@326
    31
deba@326
    32
///\ingroup matching
deba@326
    33
///\file
deba@327
    34
///\brief Maximum matching algorithms in general graphs.
deba@326
    35
deba@326
    36
namespace lemon {
deba@326
    37
deba@327
    38
  /// \ingroup matching
deba@326
    39
  ///
deba@327
    40
  /// \brief Edmonds' alternating forest maximum matching algorithm.
deba@326
    41
  ///
alpar@330
    42
  /// This class implements Edmonds' alternating forest matching
alpar@330
    43
  /// algorithm. The algorithm can be started from an arbitrary initial
alpar@330
    44
  /// matching (the default is the empty one)
deba@326
    45
  ///
alpar@330
    46
  /// The dual solution of the problem is a map of the nodes to
deba@327
    47
  /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
deba@327
    48
  /// MATCHED/C showing the Gallai-Edmonds decomposition of the
deba@327
    49
  /// graph. The nodes in \c EVEN/D induce a graph with
deba@327
    50
  /// factor-critical components, the nodes in \c ODD/A form the
deba@327
    51
  /// barrier, and the nodes in \c MATCHED/C induce a graph having a
alpar@330
    52
  /// perfect matching. The number of the factor-critical components
deba@327
    53
  /// minus the number of barrier nodes is a lower bound on the
alpar@330
    54
  /// unmatched nodes, and the matching is optimal if and only if this bound is
deba@327
    55
  /// tight. This decomposition can be attained by calling \c
deba@327
    56
  /// decomposition() after running the algorithm.
deba@326
    57
  ///
kpeter@559
    58
  /// \param GR The graph type the algorithm runs on.
kpeter@559
    59
  template <typename GR>
deba@326
    60
  class MaxMatching {
deba@327
    61
  public:
deba@327
    62
kpeter@559
    63
    typedef GR Graph;
deba@327
    64
    typedef typename Graph::template NodeMap<typename Graph::Arc>
deba@327
    65
    MatchingMap;
deba@327
    66
deba@327
    67
    ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
deba@327
    68
    ///
alpar@330
    69
    ///Indicates the Gallai-Edmonds decomposition of the graph. The
deba@327
    70
    ///nodes with Status \c EVEN/D induce a graph with factor-critical
deba@327
    71
    ///components, the nodes in \c ODD/A form the canonical barrier,
deba@327
    72
    ///and the nodes in \c MATCHED/C induce a graph having a perfect
deba@327
    73
    ///matching.
deba@327
    74
    enum Status {
deba@327
    75
      EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
deba@327
    76
    };
deba@327
    77
deba@327
    78
    typedef typename Graph::template NodeMap<Status> StatusMap;
deba@327
    79
deba@327
    80
  private:
deba@326
    81
deba@326
    82
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
deba@326
    83
deba@327
    84
    typedef UnionFindEnum<IntNodeMap> BlossomSet;
deba@327
    85
    typedef ExtendFindEnum<IntNodeMap> TreeSet;
deba@327
    86
    typedef RangeMap<Node> NodeIntMap;
deba@327
    87
    typedef MatchingMap EarMap;
deba@327
    88
    typedef std::vector<Node> NodeQueue;
deba@327
    89
deba@327
    90
    const Graph& _graph;
deba@327
    91
    MatchingMap* _matching;
deba@327
    92
    StatusMap* _status;
deba@327
    93
deba@327
    94
    EarMap* _ear;
deba@327
    95
deba@327
    96
    IntNodeMap* _blossom_set_index;
deba@327
    97
    BlossomSet* _blossom_set;
deba@327
    98
    NodeIntMap* _blossom_rep;
deba@327
    99
deba@327
   100
    IntNodeMap* _tree_set_index;
deba@327
   101
    TreeSet* _tree_set;
deba@327
   102
deba@327
   103
    NodeQueue _node_queue;
deba@327
   104
    int _process, _postpone, _last;
deba@327
   105
deba@327
   106
    int _node_num;
deba@327
   107
deba@327
   108
  private:
deba@327
   109
deba@327
   110
    void createStructures() {
deba@327
   111
      _node_num = countNodes(_graph);
deba@327
   112
      if (!_matching) {
deba@327
   113
        _matching = new MatchingMap(_graph);
deba@327
   114
      }
deba@327
   115
      if (!_status) {
deba@327
   116
        _status = new StatusMap(_graph);
deba@327
   117
      }
deba@327
   118
      if (!_ear) {
deba@327
   119
        _ear = new EarMap(_graph);
deba@327
   120
      }
deba@327
   121
      if (!_blossom_set) {
deba@327
   122
        _blossom_set_index = new IntNodeMap(_graph);
deba@327
   123
        _blossom_set = new BlossomSet(*_blossom_set_index);
deba@327
   124
      }
deba@327
   125
      if (!_blossom_rep) {
deba@327
   126
        _blossom_rep = new NodeIntMap(_node_num);
deba@327
   127
      }
deba@327
   128
      if (!_tree_set) {
deba@327
   129
        _tree_set_index = new IntNodeMap(_graph);
deba@327
   130
        _tree_set = new TreeSet(*_tree_set_index);
deba@327
   131
      }
deba@327
   132
      _node_queue.resize(_node_num);
deba@327
   133
    }
deba@327
   134
deba@327
   135
    void destroyStructures() {
deba@327
   136
      if (_matching) {
deba@327
   137
        delete _matching;
deba@327
   138
      }
deba@327
   139
      if (_status) {
deba@327
   140
        delete _status;
deba@327
   141
      }
deba@327
   142
      if (_ear) {
deba@327
   143
        delete _ear;
deba@327
   144
      }
deba@327
   145
      if (_blossom_set) {
deba@327
   146
        delete _blossom_set;
deba@327
   147
        delete _blossom_set_index;
deba@327
   148
      }
deba@327
   149
      if (_blossom_rep) {
deba@327
   150
        delete _blossom_rep;
deba@327
   151
      }
deba@327
   152
      if (_tree_set) {
deba@327
   153
        delete _tree_set_index;
deba@327
   154
        delete _tree_set;
deba@327
   155
      }
deba@327
   156
    }
deba@327
   157
deba@327
   158
    void processDense(const Node& n) {
deba@327
   159
      _process = _postpone = _last = 0;
deba@327
   160
      _node_queue[_last++] = n;
deba@327
   161
deba@327
   162
      while (_process != _last) {
deba@327
   163
        Node u = _node_queue[_process++];
deba@327
   164
        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
deba@327
   165
          Node v = _graph.target(a);
deba@327
   166
          if ((*_status)[v] == MATCHED) {
deba@327
   167
            extendOnArc(a);
deba@327
   168
          } else if ((*_status)[v] == UNMATCHED) {
deba@327
   169
            augmentOnArc(a);
deba@327
   170
            return;
deba@327
   171
          }
deba@327
   172
        }
deba@327
   173
      }
deba@327
   174
deba@327
   175
      while (_postpone != _last) {
deba@327
   176
        Node u = _node_queue[_postpone++];
deba@327
   177
deba@327
   178
        for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
deba@327
   179
          Node v = _graph.target(a);
deba@327
   180
deba@327
   181
          if ((*_status)[v] == EVEN) {
deba@327
   182
            if (_blossom_set->find(u) != _blossom_set->find(v)) {
deba@327
   183
              shrinkOnEdge(a);
deba@327
   184
            }
deba@327
   185
          }
deba@327
   186
deba@327
   187
          while (_process != _last) {
deba@327
   188
            Node w = _node_queue[_process++];
deba@327
   189
            for (OutArcIt b(_graph, w); b != INVALID; ++b) {
deba@327
   190
              Node x = _graph.target(b);
deba@327
   191
              if ((*_status)[x] == MATCHED) {
deba@327
   192
                extendOnArc(b);
deba@327
   193
              } else if ((*_status)[x] == UNMATCHED) {
deba@327
   194
                augmentOnArc(b);
deba@327
   195
                return;
deba@327
   196
              }
deba@327
   197
            }
deba@327
   198
          }
deba@327
   199
        }
deba@327
   200
      }
deba@327
   201
    }
deba@327
   202
deba@327
   203
    void processSparse(const Node& n) {
deba@327
   204
      _process = _last = 0;
deba@327
   205
      _node_queue[_last++] = n;
deba@327
   206
      while (_process != _last) {
deba@327
   207
        Node u = _node_queue[_process++];
deba@327
   208
        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
deba@327
   209
          Node v = _graph.target(a);
deba@327
   210
deba@327
   211
          if ((*_status)[v] == EVEN) {
deba@327
   212
            if (_blossom_set->find(u) != _blossom_set->find(v)) {
deba@327
   213
              shrinkOnEdge(a);
deba@327
   214
            }
deba@327
   215
          } else if ((*_status)[v] == MATCHED) {
deba@327
   216
            extendOnArc(a);
deba@327
   217
          } else if ((*_status)[v] == UNMATCHED) {
deba@327
   218
            augmentOnArc(a);
deba@327
   219
            return;
deba@327
   220
          }
deba@327
   221
        }
deba@327
   222
      }
deba@327
   223
    }
deba@327
   224
deba@327
   225
    void shrinkOnEdge(const Edge& e) {
deba@327
   226
      Node nca = INVALID;
deba@327
   227
deba@327
   228
      {
deba@327
   229
        std::set<Node> left_set, right_set;
deba@327
   230
deba@327
   231
        Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
deba@327
   232
        left_set.insert(left);
deba@327
   233
deba@327
   234
        Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
deba@327
   235
        right_set.insert(right);
deba@327
   236
deba@327
   237
        while (true) {
deba@327
   238
          if ((*_matching)[left] == INVALID) break;
deba@327
   239
          left = _graph.target((*_matching)[left]);
deba@327
   240
          left = (*_blossom_rep)[_blossom_set->
deba@327
   241
                                 find(_graph.target((*_ear)[left]))];
deba@327
   242
          if (right_set.find(left) != right_set.end()) {
deba@327
   243
            nca = left;
deba@327
   244
            break;
deba@327
   245
          }
deba@327
   246
          left_set.insert(left);
deba@327
   247
deba@327
   248
          if ((*_matching)[right] == INVALID) break;
deba@327
   249
          right = _graph.target((*_matching)[right]);
deba@327
   250
          right = (*_blossom_rep)[_blossom_set->
deba@327
   251
                                  find(_graph.target((*_ear)[right]))];
deba@327
   252
          if (left_set.find(right) != left_set.end()) {
deba@327
   253
            nca = right;
deba@327
   254
            break;
deba@327
   255
          }
deba@327
   256
          right_set.insert(right);
deba@327
   257
        }
deba@327
   258
deba@327
   259
        if (nca == INVALID) {
deba@327
   260
          if ((*_matching)[left] == INVALID) {
deba@327
   261
            nca = right;
deba@327
   262
            while (left_set.find(nca) == left_set.end()) {
deba@327
   263
              nca = _graph.target((*_matching)[nca]);
deba@327
   264
              nca =(*_blossom_rep)[_blossom_set->
deba@327
   265
                                   find(_graph.target((*_ear)[nca]))];
deba@327
   266
            }
deba@327
   267
          } else {
deba@327
   268
            nca = left;
deba@327
   269
            while (right_set.find(nca) == right_set.end()) {
deba@327
   270
              nca = _graph.target((*_matching)[nca]);
deba@327
   271
              nca = (*_blossom_rep)[_blossom_set->
deba@327
   272
                                   find(_graph.target((*_ear)[nca]))];
deba@327
   273
            }
deba@327
   274
          }
deba@327
   275
        }
deba@327
   276
      }
deba@327
   277
deba@327
   278
      {
deba@327
   279
deba@327
   280
        Node node = _graph.u(e);
deba@327
   281
        Arc arc = _graph.direct(e, true);
deba@327
   282
        Node base = (*_blossom_rep)[_blossom_set->find(node)];
deba@327
   283
deba@327
   284
        while (base != nca) {
kpeter@581
   285
          (*_ear)[node] = arc;
deba@327
   286
deba@327
   287
          Node n = node;
deba@327
   288
          while (n != base) {
deba@327
   289
            n = _graph.target((*_matching)[n]);
deba@327
   290
            Arc a = (*_ear)[n];
deba@327
   291
            n = _graph.target(a);
kpeter@581
   292
            (*_ear)[n] = _graph.oppositeArc(a);
deba@327
   293
          }
deba@327
   294
          node = _graph.target((*_matching)[base]);
deba@327
   295
          _tree_set->erase(base);
deba@327
   296
          _tree_set->erase(node);
deba@327
   297
          _blossom_set->insert(node, _blossom_set->find(base));
kpeter@581
   298
          (*_status)[node] = EVEN;
deba@327
   299
          _node_queue[_last++] = node;
deba@327
   300
          arc = _graph.oppositeArc((*_ear)[node]);
deba@327
   301
          node = _graph.target((*_ear)[node]);
deba@327
   302
          base = (*_blossom_rep)[_blossom_set->find(node)];
deba@327
   303
          _blossom_set->join(_graph.target(arc), base);
deba@327
   304
        }
deba@327
   305
      }
deba@327
   306
kpeter@581
   307
      (*_blossom_rep)[_blossom_set->find(nca)] = nca;
deba@327
   308
deba@327
   309
      {
deba@327
   310
deba@327
   311
        Node node = _graph.v(e);
deba@327
   312
        Arc arc = _graph.direct(e, false);
deba@327
   313
        Node base = (*_blossom_rep)[_blossom_set->find(node)];
deba@327
   314
deba@327
   315
        while (base != nca) {
kpeter@581
   316
          (*_ear)[node] = arc;
deba@327
   317
deba@327
   318
          Node n = node;
deba@327
   319
          while (n != base) {
deba@327
   320
            n = _graph.target((*_matching)[n]);
deba@327
   321
            Arc a = (*_ear)[n];
deba@327
   322
            n = _graph.target(a);
kpeter@581
   323
            (*_ear)[n] = _graph.oppositeArc(a);
deba@327
   324
          }
deba@327
   325
          node = _graph.target((*_matching)[base]);
deba@327
   326
          _tree_set->erase(base);
deba@327
   327
          _tree_set->erase(node);
deba@327
   328
          _blossom_set->insert(node, _blossom_set->find(base));
kpeter@581
   329
          (*_status)[node] = EVEN;
deba@327
   330
          _node_queue[_last++] = node;
deba@327
   331
          arc = _graph.oppositeArc((*_ear)[node]);
deba@327
   332
          node = _graph.target((*_ear)[node]);
deba@327
   333
          base = (*_blossom_rep)[_blossom_set->find(node)];
deba@327
   334
          _blossom_set->join(_graph.target(arc), base);
deba@327
   335
        }
deba@327
   336
      }
deba@327
   337
kpeter@581
   338
      (*_blossom_rep)[_blossom_set->find(nca)] = nca;
deba@327
   339
    }
deba@327
   340
deba@327
   341
deba@327
   342
deba@327
   343
    void extendOnArc(const Arc& a) {
deba@327
   344
      Node base = _graph.source(a);
deba@327
   345
      Node odd = _graph.target(a);
deba@327
   346
kpeter@581
   347
      (*_ear)[odd] = _graph.oppositeArc(a);
deba@327
   348
      Node even = _graph.target((*_matching)[odd]);
kpeter@581
   349
      (*_blossom_rep)[_blossom_set->insert(even)] = even;
kpeter@581
   350
      (*_status)[odd] = ODD;
kpeter@581
   351
      (*_status)[even] = EVEN;
deba@327
   352
      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
deba@327
   353
      _tree_set->insert(odd, tree);
deba@327
   354
      _tree_set->insert(even, tree);
deba@327
   355
      _node_queue[_last++] = even;
deba@327
   356
deba@327
   357
    }
deba@327
   358
deba@327
   359
    void augmentOnArc(const Arc& a) {
deba@327
   360
      Node even = _graph.source(a);
deba@327
   361
      Node odd = _graph.target(a);
deba@327
   362
deba@327
   363
      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
deba@327
   364
kpeter@581
   365
      (*_matching)[odd] = _graph.oppositeArc(a);
kpeter@581
   366
      (*_status)[odd] = MATCHED;
deba@327
   367
deba@327
   368
      Arc arc = (*_matching)[even];
kpeter@581
   369
      (*_matching)[even] = a;
deba@327
   370
deba@327
   371
      while (arc != INVALID) {
deba@327
   372
        odd = _graph.target(arc);
deba@327
   373
        arc = (*_ear)[odd];
deba@327
   374
        even = _graph.target(arc);
kpeter@581
   375
        (*_matching)[odd] = arc;
deba@327
   376
        arc = (*_matching)[even];
kpeter@581
   377
        (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]);
deba@327
   378
      }
deba@327
   379
deba@327
   380
      for (typename TreeSet::ItemIt it(*_tree_set, tree);
deba@327
   381
           it != INVALID; ++it) {
deba@327
   382
        if ((*_status)[it] == ODD) {
kpeter@581
   383
          (*_status)[it] = MATCHED;
deba@327
   384
        } else {
deba@327
   385
          int blossom = _blossom_set->find(it);
deba@327
   386
          for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
deba@327
   387
               jt != INVALID; ++jt) {
kpeter@581
   388
            (*_status)[jt] = MATCHED;
deba@327
   389
          }
deba@327
   390
          _blossom_set->eraseClass(blossom);
deba@327
   391
        }
deba@327
   392
      }
deba@327
   393
      _tree_set->eraseClass(tree);
deba@327
   394
deba@327
   395
    }
deba@326
   396
deba@326
   397
  public:
deba@326
   398
deba@327
   399
    /// \brief Constructor
deba@326
   400
    ///
deba@327
   401
    /// Constructor.
deba@327
   402
    MaxMatching(const Graph& graph)
deba@327
   403
      : _graph(graph), _matching(0), _status(0), _ear(0),
deba@327
   404
        _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
deba@327
   405
        _tree_set_index(0), _tree_set(0) {}
deba@327
   406
deba@327
   407
    ~MaxMatching() {
deba@327
   408
      destroyStructures();
deba@327
   409
    }
deba@327
   410
deba@327
   411
    /// \name Execution control
alpar@330
   412
    /// The simplest way to execute the algorithm is to use the
deba@327
   413
    /// \c run() member function.
deba@327
   414
    /// \n
deba@327
   415
alpar@330
   416
    /// If you need better control on the execution, you must call
deba@327
   417
    /// \ref init(), \ref greedyInit() or \ref matchingInit()
alpar@330
   418
    /// functions first, then you can start the algorithm with the \ref
deba@425
   419
    /// startSparse() or startDense() functions.
deba@327
   420
deba@327
   421
    ///@{
deba@327
   422
deba@327
   423
    /// \brief Sets the actual matching to the empty matching.
deba@326
   424
    ///
deba@327
   425
    /// Sets the actual matching to the empty matching.
deba@326
   426
    ///
deba@326
   427
    void init() {
deba@327
   428
      createStructures();
deba@327
   429
      for(NodeIt n(_graph); n != INVALID; ++n) {
kpeter@581
   430
        (*_matching)[n] = INVALID;
kpeter@581
   431
        (*_status)[n] = UNMATCHED;
deba@326
   432
      }
deba@326
   433
    }
deba@326
   434
alpar@330
   435
    ///\brief Finds an initial matching in a greedy way
deba@326
   436
    ///
alpar@330
   437
    ///It finds an initial matching in a greedy way.
deba@326
   438
    void greedyInit() {
deba@327
   439
      createStructures();
deba@327
   440
      for (NodeIt n(_graph); n != INVALID; ++n) {
kpeter@581
   441
        (*_matching)[n] = INVALID;
kpeter@581
   442
        (*_status)[n] = UNMATCHED;
deba@326
   443
      }
deba@327
   444
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@327
   445
        if ((*_matching)[n] == INVALID) {
deba@327
   446
          for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
deba@327
   447
            Node v = _graph.target(a);
deba@327
   448
            if ((*_matching)[v] == INVALID && v != n) {
kpeter@581
   449
              (*_matching)[n] = a;
kpeter@581
   450
              (*_status)[n] = MATCHED;
kpeter@581
   451
              (*_matching)[v] = _graph.oppositeArc(a);
kpeter@581
   452
              (*_status)[v] = MATCHED;
deba@326
   453
              break;
deba@326
   454
            }
deba@326
   455
          }
deba@326
   456
        }
deba@326
   457
      }
deba@326
   458
    }
deba@326
   459
deba@327
   460
alpar@330
   461
    /// \brief Initialize the matching from a map containing.
deba@326
   462
    ///
deba@327
   463
    /// Initialize the matching from a \c bool valued \c Edge map. This
deba@327
   464
    /// map must have the property that there are no two incident edges
deba@327
   465
    /// with true value, ie. it contains a matching.
kpeter@559
   466
    /// \return \c true if the map contains a matching.
deba@327
   467
    template <typename MatchingMap>
deba@327
   468
    bool matchingInit(const MatchingMap& matching) {
deba@327
   469
      createStructures();
deba@327
   470
deba@327
   471
      for (NodeIt n(_graph); n != INVALID; ++n) {
kpeter@581
   472
        (*_matching)[n] = INVALID;
kpeter@581
   473
        (*_status)[n] = UNMATCHED;
deba@326
   474
      }
deba@327
   475
      for(EdgeIt e(_graph); e!=INVALID; ++e) {
deba@327
   476
        if (matching[e]) {
deba@327
   477
deba@327
   478
          Node u = _graph.u(e);
deba@327
   479
          if ((*_matching)[u] != INVALID) return false;
kpeter@581
   480
          (*_matching)[u] = _graph.direct(e, true);
kpeter@581
   481
          (*_status)[u] = MATCHED;
deba@327
   482
deba@327
   483
          Node v = _graph.v(e);
deba@327
   484
          if ((*_matching)[v] != INVALID) return false;
kpeter@581
   485
          (*_matching)[v] = _graph.direct(e, false);
kpeter@581
   486
          (*_status)[v] = MATCHED;
deba@327
   487
        }
deba@327
   488
      }
deba@327
   489
      return true;
deba@326
   490
    }
deba@326
   491
deba@327
   492
    /// \brief Starts Edmonds' algorithm
deba@326
   493
    ///
deba@327
   494
    /// If runs the original Edmonds' algorithm.
deba@327
   495
    void startSparse() {
deba@327
   496
      for(NodeIt n(_graph); n != INVALID; ++n) {
deba@327
   497
        if ((*_status)[n] == UNMATCHED) {
deba@327
   498
          (*_blossom_rep)[_blossom_set->insert(n)] = n;
deba@327
   499
          _tree_set->insert(n);
kpeter@581
   500
          (*_status)[n] = EVEN;
deba@327
   501
          processSparse(n);
deba@326
   502
        }
deba@326
   503
      }
deba@326
   504
    }
deba@326
   505
deba@327
   506
    /// \brief Starts Edmonds' algorithm.
deba@326
   507
    ///
deba@327
   508
    /// It runs Edmonds' algorithm with a heuristic of postponing
alpar@330
   509
    /// shrinks, therefore resulting in a faster algorithm for dense graphs.
deba@327
   510
    void startDense() {
deba@327
   511
      for(NodeIt n(_graph); n != INVALID; ++n) {
deba@327
   512
        if ((*_status)[n] == UNMATCHED) {
deba@327
   513
          (*_blossom_rep)[_blossom_set->insert(n)] = n;
deba@327
   514
          _tree_set->insert(n);
kpeter@581
   515
          (*_status)[n] = EVEN;
deba@327
   516
          processDense(n);
deba@327
   517
        }
deba@327
   518
      }
deba@327
   519
    }
deba@327
   520
deba@327
   521
deba@327
   522
    /// \brief Runs Edmonds' algorithm
deba@327
   523
    ///
deba@327
   524
    /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
deba@327
   525
    /// or Edmonds' algorithm with a heuristic of
deba@327
   526
    /// postponing shrinks for dense graphs.
deba@326
   527
    void run() {
deba@327
   528
      if (countEdges(_graph) < 2 * countNodes(_graph)) {
deba@326
   529
        greedyInit();
deba@326
   530
        startSparse();
deba@326
   531
      } else {
deba@326
   532
        init();
deba@326
   533
        startDense();
deba@326
   534
      }
deba@326
   535
    }
deba@326
   536
deba@327
   537
    /// @}
deba@327
   538
deba@327
   539
    /// \name Primal solution
alpar@330
   540
    /// Functions to get the primal solution, ie. the matching.
deba@327
   541
deba@327
   542
    /// @{
deba@326
   543
alpar@330
   544
    ///\brief Returns the size of the current matching.
deba@326
   545
    ///
alpar@330
   546
    ///Returns the size of the current matching. After \ref
deba@327
   547
    ///run() it returns the size of the maximum matching in the graph.
deba@327
   548
    int matchingSize() const {
deba@327
   549
      int size = 0;
deba@327
   550
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@327
   551
        if ((*_matching)[n] != INVALID) {
deba@327
   552
          ++size;
deba@326
   553
        }
deba@326
   554
      }
deba@327
   555
      return size / 2;
deba@326
   556
    }
deba@326
   557
deba@327
   558
    /// \brief Returns true when the edge is in the matching.
deba@327
   559
    ///
deba@327
   560
    /// Returns true when the edge is in the matching.
deba@327
   561
    bool matching(const Edge& edge) const {
deba@327
   562
      return edge == (*_matching)[_graph.u(edge)];
deba@327
   563
    }
deba@327
   564
deba@327
   565
    /// \brief Returns the matching edge incident to the given node.
deba@327
   566
    ///
deba@327
   567
    /// Returns the matching edge of a \c node in the actual matching or
deba@327
   568
    /// INVALID if the \c node is not covered by the actual matching.
deba@327
   569
    Arc matching(const Node& n) const {
deba@327
   570
      return (*_matching)[n];
deba@327
   571
    }
deba@326
   572
deba@326
   573
    ///\brief Returns the mate of a node in the actual matching.
deba@326
   574
    ///
deba@327
   575
    ///Returns the mate of a \c node in the actual matching or
deba@327
   576
    ///INVALID if the \c node is not covered by the actual matching.
deba@327
   577
    Node mate(const Node& n) const {
deba@327
   578
      return (*_matching)[n] != INVALID ?
deba@327
   579
        _graph.target((*_matching)[n]) : INVALID;
deba@326
   580
    }
deba@326
   581
deba@327
   582
    /// @}
deba@327
   583
deba@327
   584
    /// \name Dual solution
alpar@330
   585
    /// Functions to get the dual solution, ie. the decomposition.
deba@327
   586
deba@327
   587
    /// @{
deba@326
   588
deba@326
   589
    /// \brief Returns the class of the node in the Edmonds-Gallai
deba@326
   590
    /// decomposition.
deba@326
   591
    ///
deba@326
   592
    /// Returns the class of the node in the Edmonds-Gallai
deba@326
   593
    /// decomposition.
deba@327
   594
    Status decomposition(const Node& n) const {
deba@327
   595
      return (*_status)[n];
deba@326
   596
    }
deba@326
   597
deba@326
   598
    /// \brief Returns true when the node is in the barrier.
deba@326
   599
    ///
deba@326
   600
    /// Returns true when the node is in the barrier.
deba@327
   601
    bool barrier(const Node& n) const {
deba@327
   602
      return (*_status)[n] == ODD;
deba@326
   603
    }
deba@326
   604
deba@327
   605
    /// @}
deba@326
   606
deba@326
   607
  };
deba@326
   608
deba@326
   609
  /// \ingroup matching
deba@326
   610
  ///
deba@326
   611
  /// \brief Weighted matching in general graphs
deba@326
   612
  ///
deba@326
   613
  /// This class provides an efficient implementation of Edmond's
deba@326
   614
  /// maximum weighted matching algorithm. The implementation is based
deba@326
   615
  /// on extensive use of priority queues and provides
kpeter@559
   616
  /// \f$O(nm\log n)\f$ time complexity.
deba@326
   617
  ///
deba@326
   618
  /// The maximum weighted matching problem is to find undirected
deba@327
   619
  /// edges in the graph with maximum overall weight and no two of
deba@327
   620
  /// them shares their ends. The problem can be formulated with the
deba@327
   621
  /// following linear program.
deba@326
   622
  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
deba@327
   623
  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
deba@327
   624
      \quad \forall B\in\mathcal{O}\f] */
deba@326
   625
  /// \f[x_e \ge 0\quad \forall e\in E\f]
deba@326
   626
  /// \f[\max \sum_{e\in E}x_ew_e\f]
deba@327
   627
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
deba@327
   628
  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
deba@327
   629
  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
deba@327
   630
  /// subsets of the nodes.
deba@326
   631
  ///
deba@326
   632
  /// The algorithm calculates an optimal matching and a proof of the
deba@326
   633
  /// optimality. The solution of the dual problem can be used to check
deba@327
   634
  /// the result of the algorithm. The dual linear problem is the
deba@327
   635
  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
deba@327
   636
      z_B \ge w_{uv} \quad \forall uv\in E\f] */
deba@326
   637
  /// \f[y_u \ge 0 \quad \forall u \in V\f]
deba@326
   638
  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
deba@327
   639
  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
deba@327
   640
      \frac{\vert B \vert - 1}{2}z_B\f] */
deba@326
   641
  ///
deba@326
   642
  /// The algorithm can be executed with \c run() or the \c init() and
deba@326
   643
  /// then the \c start() member functions. After it the matching can
deba@326
   644
  /// be asked with \c matching() or mate() functions. The dual
deba@326
   645
  /// solution can be get with \c nodeValue(), \c blossomNum() and \c
deba@326
   646
  /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
alpar@330
   647
  /// "BlossomIt" nested class, which is able to iterate on the nodes
deba@326
   648
  /// of a blossom. If the value type is integral then the dual
deba@326
   649
  /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
kpeter@559
   650
  template <typename GR,
kpeter@559
   651
            typename WM = typename GR::template EdgeMap<int> >
deba@326
   652
  class MaxWeightedMatching {
deba@326
   653
  public:
deba@326
   654
kpeter@559
   655
    ///\e
kpeter@559
   656
    typedef GR Graph;
kpeter@559
   657
    ///\e
kpeter@559
   658
    typedef WM WeightMap;
kpeter@559
   659
    ///\e
deba@326
   660
    typedef typename WeightMap::Value Value;
deba@326
   661
deba@326
   662
    /// \brief Scaling factor for dual solution
deba@326
   663
    ///
deba@326
   664
    /// Scaling factor for dual solution, it is equal to 4 or 1
deba@326
   665
    /// according to the value type.
deba@326
   666
    static const int dualScale =
deba@326
   667
      std::numeric_limits<Value>::is_integer ? 4 : 1;
deba@326
   668
deba@326
   669
    typedef typename Graph::template NodeMap<typename Graph::Arc>
deba@326
   670
    MatchingMap;
deba@326
   671
deba@326
   672
  private:
deba@326
   673
deba@326
   674
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
deba@326
   675
deba@326
   676
    typedef typename Graph::template NodeMap<Value> NodePotential;
deba@326
   677
    typedef std::vector<Node> BlossomNodeList;
deba@326
   678
deba@326
   679
    struct BlossomVariable {
deba@326
   680
      int begin, end;
deba@326
   681
      Value value;
deba@326
   682
deba@326
   683
      BlossomVariable(int _begin, int _end, Value _value)
deba@326
   684
        : begin(_begin), end(_end), value(_value) {}
deba@326
   685
deba@326
   686
    };
deba@326
   687
deba@326
   688
    typedef std::vector<BlossomVariable> BlossomPotential;
deba@326
   689
deba@326
   690
    const Graph& _graph;
deba@326
   691
    const WeightMap& _weight;
deba@326
   692
deba@326
   693
    MatchingMap* _matching;
deba@326
   694
deba@326
   695
    NodePotential* _node_potential;
deba@326
   696
deba@326
   697
    BlossomPotential _blossom_potential;
deba@326
   698
    BlossomNodeList _blossom_node_list;
deba@326
   699
deba@326
   700
    int _node_num;
deba@326
   701
    int _blossom_num;
deba@326
   702
deba@326
   703
    typedef RangeMap<int> IntIntMap;
deba@326
   704
deba@326
   705
    enum Status {
deba@326
   706
      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
deba@326
   707
    };
deba@326
   708
deba@327
   709
    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
deba@326
   710
    struct BlossomData {
deba@326
   711
      int tree;
deba@326
   712
      Status status;
deba@326
   713
      Arc pred, next;
deba@326
   714
      Value pot, offset;
deba@326
   715
      Node base;
deba@326
   716
    };
deba@326
   717
deba@327
   718
    IntNodeMap *_blossom_index;
deba@326
   719
    BlossomSet *_blossom_set;
deba@326
   720
    RangeMap<BlossomData>* _blossom_data;
deba@326
   721
deba@327
   722
    IntNodeMap *_node_index;
deba@327
   723
    IntArcMap *_node_heap_index;
deba@326
   724
deba@326
   725
    struct NodeData {
deba@326
   726
deba@327
   727
      NodeData(IntArcMap& node_heap_index)
deba@326
   728
        : heap(node_heap_index) {}
deba@326
   729
deba@326
   730
      int blossom;
deba@326
   731
      Value pot;
deba@327
   732
      BinHeap<Value, IntArcMap> heap;
deba@326
   733
      std::map<int, Arc> heap_index;
deba@326
   734
deba@326
   735
      int tree;
deba@326
   736
    };
deba@326
   737
deba@326
   738
    RangeMap<NodeData>* _node_data;
deba@326
   739
deba@326
   740
    typedef ExtendFindEnum<IntIntMap> TreeSet;
deba@326
   741
deba@326
   742
    IntIntMap *_tree_set_index;
deba@326
   743
    TreeSet *_tree_set;
deba@326
   744
deba@327
   745
    IntNodeMap *_delta1_index;
deba@327
   746
    BinHeap<Value, IntNodeMap> *_delta1;
deba@326
   747
deba@326
   748
    IntIntMap *_delta2_index;
deba@326
   749
    BinHeap<Value, IntIntMap> *_delta2;
deba@326
   750
deba@327
   751
    IntEdgeMap *_delta3_index;
deba@327
   752
    BinHeap<Value, IntEdgeMap> *_delta3;
deba@326
   753
deba@326
   754
    IntIntMap *_delta4_index;
deba@326
   755
    BinHeap<Value, IntIntMap> *_delta4;
deba@326
   756
deba@326
   757
    Value _delta_sum;
deba@326
   758
deba@326
   759
    void createStructures() {
deba@326
   760
      _node_num = countNodes(_graph);
deba@326
   761
      _blossom_num = _node_num * 3 / 2;
deba@326
   762
deba@326
   763
      if (!_matching) {
deba@326
   764
        _matching = new MatchingMap(_graph);
deba@326
   765
      }
deba@326
   766
      if (!_node_potential) {
deba@326
   767
        _node_potential = new NodePotential(_graph);
deba@326
   768
      }
deba@326
   769
      if (!_blossom_set) {
deba@327
   770
        _blossom_index = new IntNodeMap(_graph);
deba@326
   771
        _blossom_set = new BlossomSet(*_blossom_index);
deba@326
   772
        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
deba@326
   773
      }
deba@326
   774
deba@326
   775
      if (!_node_index) {
deba@327
   776
        _node_index = new IntNodeMap(_graph);
deba@327
   777
        _node_heap_index = new IntArcMap(_graph);
deba@326
   778
        _node_data = new RangeMap<NodeData>(_node_num,
deba@326
   779
                                              NodeData(*_node_heap_index));
deba@326
   780
      }
deba@326
   781
deba@326
   782
      if (!_tree_set) {
deba@326
   783
        _tree_set_index = new IntIntMap(_blossom_num);
deba@326
   784
        _tree_set = new TreeSet(*_tree_set_index);
deba@326
   785
      }
deba@326
   786
      if (!_delta1) {
deba@327
   787
        _delta1_index = new IntNodeMap(_graph);
deba@327
   788
        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
deba@326
   789
      }
deba@326
   790
      if (!_delta2) {
deba@326
   791
        _delta2_index = new IntIntMap(_blossom_num);
deba@326
   792
        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
deba@326
   793
      }
deba@326
   794
      if (!_delta3) {
deba@327
   795
        _delta3_index = new IntEdgeMap(_graph);
deba@327
   796
        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
deba@326
   797
      }
deba@326
   798
      if (!_delta4) {
deba@326
   799
        _delta4_index = new IntIntMap(_blossom_num);
deba@326
   800
        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
deba@326
   801
      }
deba@326
   802
    }
deba@326
   803
deba@326
   804
    void destroyStructures() {
deba@326
   805
      _node_num = countNodes(_graph);
deba@326
   806
      _blossom_num = _node_num * 3 / 2;
deba@326
   807
deba@326
   808
      if (_matching) {
deba@326
   809
        delete _matching;
deba@326
   810
      }
deba@326
   811
      if (_node_potential) {
deba@326
   812
        delete _node_potential;
deba@326
   813
      }
deba@326
   814
      if (_blossom_set) {
deba@326
   815
        delete _blossom_index;
deba@326
   816
        delete _blossom_set;
deba@326
   817
        delete _blossom_data;
deba@326
   818
      }
deba@326
   819
deba@326
   820
      if (_node_index) {
deba@326
   821
        delete _node_index;
deba@326
   822
        delete _node_heap_index;
deba@326
   823
        delete _node_data;
deba@326
   824
      }
deba@326
   825
deba@326
   826
      if (_tree_set) {
deba@326
   827
        delete _tree_set_index;
deba@326
   828
        delete _tree_set;
deba@326
   829
      }
deba@326
   830
      if (_delta1) {
deba@326
   831
        delete _delta1_index;
deba@326
   832
        delete _delta1;
deba@326
   833
      }
deba@326
   834
      if (_delta2) {
deba@326
   835
        delete _delta2_index;
deba@326
   836
        delete _delta2;
deba@326
   837
      }
deba@326
   838
      if (_delta3) {
deba@326
   839
        delete _delta3_index;
deba@326
   840
        delete _delta3;
deba@326
   841
      }
deba@326
   842
      if (_delta4) {
deba@326
   843
        delete _delta4_index;
deba@326
   844
        delete _delta4;
deba@326
   845
      }
deba@326
   846
    }
deba@326
   847
deba@326
   848
    void matchedToEven(int blossom, int tree) {
deba@326
   849
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
deba@326
   850
        _delta2->erase(blossom);
deba@326
   851
      }
deba@326
   852
deba@326
   853
      if (!_blossom_set->trivial(blossom)) {
deba@326
   854
        (*_blossom_data)[blossom].pot -=
deba@326
   855
          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
deba@326
   856
      }
deba@326
   857
deba@326
   858
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
deba@326
   859
           n != INVALID; ++n) {
deba@326
   860
deba@326
   861
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
deba@326
   862
        int ni = (*_node_index)[n];
deba@326
   863
deba@326
   864
        (*_node_data)[ni].heap.clear();
deba@326
   865
        (*_node_data)[ni].heap_index.clear();
deba@326
   866
deba@326
   867
        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
deba@326
   868
deba@326
   869
        _delta1->push(n, (*_node_data)[ni].pot);
deba@326
   870
deba@326
   871
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
deba@326
   872
          Node v = _graph.source(e);
deba@326
   873
          int vb = _blossom_set->find(v);
deba@326
   874
          int vi = (*_node_index)[v];
deba@326
   875
deba@326
   876
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
deba@326
   877
            dualScale * _weight[e];
deba@326
   878
deba@326
   879
          if ((*_blossom_data)[vb].status == EVEN) {
deba@326
   880
            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
deba@326
   881
              _delta3->push(e, rw / 2);
deba@326
   882
            }
deba@326
   883
          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
deba@326
   884
            if (_delta3->state(e) != _delta3->IN_HEAP) {
deba@326
   885
              _delta3->push(e, rw);
deba@326
   886
            }
deba@326
   887
          } else {
deba@326
   888
            typename std::map<int, Arc>::iterator it =
deba@326
   889
              (*_node_data)[vi].heap_index.find(tree);
deba@326
   890
deba@326
   891
            if (it != (*_node_data)[vi].heap_index.end()) {
deba@326
   892
              if ((*_node_data)[vi].heap[it->second] > rw) {
deba@326
   893
                (*_node_data)[vi].heap.replace(it->second, e);
deba@326
   894
                (*_node_data)[vi].heap.decrease(e, rw);
deba@326
   895
                it->second = e;
deba@326
   896
              }
deba@326
   897
            } else {
deba@326
   898
              (*_node_data)[vi].heap.push(e, rw);
deba@326
   899
              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
deba@326
   900
            }
deba@326
   901
deba@326
   902
            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
deba@326
   903
              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
deba@326
   904
deba@326
   905
              if ((*_blossom_data)[vb].status == MATCHED) {
deba@326
   906
                if (_delta2->state(vb) != _delta2->IN_HEAP) {
deba@326
   907
                  _delta2->push(vb, _blossom_set->classPrio(vb) -
deba@326
   908
                               (*_blossom_data)[vb].offset);
deba@326
   909
                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
deba@326
   910
                           (*_blossom_data)[vb].offset){
deba@326
   911
                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
deba@326
   912
                                   (*_blossom_data)[vb].offset);
deba@326
   913
                }
deba@326
   914
              }
deba@326
   915
            }
deba@326
   916
          }
deba@326
   917
        }
deba@326
   918
      }
deba@326
   919
      (*_blossom_data)[blossom].offset = 0;
deba@326
   920
    }
deba@326
   921
deba@326
   922
    void matchedToOdd(int blossom) {
deba@326
   923
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
deba@326
   924
        _delta2->erase(blossom);
deba@326
   925
      }
deba@326
   926
      (*_blossom_data)[blossom].offset += _delta_sum;
deba@326
   927
      if (!_blossom_set->trivial(blossom)) {
deba@326
   928
        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
deba@326
   929
                     (*_blossom_data)[blossom].offset);
deba@326
   930
      }
deba@326
   931
    }
deba@326
   932
deba@326
   933
    void evenToMatched(int blossom, int tree) {
deba@326
   934
      if (!_blossom_set->trivial(blossom)) {
deba@326
   935
        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
deba@326
   936
      }
deba@326
   937
deba@326
   938
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
deba@326
   939
           n != INVALID; ++n) {
deba@326
   940
        int ni = (*_node_index)[n];
deba@326
   941
        (*_node_data)[ni].pot -= _delta_sum;
deba@326
   942
deba@326
   943
        _delta1->erase(n);
deba@326
   944
deba@326
   945
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
deba@326
   946
          Node v = _graph.source(e);
deba@326
   947
          int vb = _blossom_set->find(v);
deba@326
   948
          int vi = (*_node_index)[v];
deba@326
   949
deba@326
   950
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
deba@326
   951
            dualScale * _weight[e];
deba@326
   952
deba@326
   953
          if (vb == blossom) {
deba@326
   954
            if (_delta3->state(e) == _delta3->IN_HEAP) {
deba@326
   955
              _delta3->erase(e);
deba@326
   956
            }
deba@326
   957
          } else if ((*_blossom_data)[vb].status == EVEN) {
deba@326
   958
deba@326
   959
            if (_delta3->state(e) == _delta3->IN_HEAP) {
deba@326
   960
              _delta3->erase(e);
deba@326
   961
            }
deba@326
   962
deba@326
   963
            int vt = _tree_set->find(vb);
deba@326
   964
deba@326
   965
            if (vt != tree) {
deba@326
   966
deba@326
   967
              Arc r = _graph.oppositeArc(e);
deba@326
   968
deba@326
   969
              typename std::map<int, Arc>::iterator it =
deba@326
   970
                (*_node_data)[ni].heap_index.find(vt);
deba@326
   971
deba@326
   972
              if (it != (*_node_data)[ni].heap_index.end()) {
deba@326
   973
                if ((*_node_data)[ni].heap[it->second] > rw) {
deba@326
   974
                  (*_node_data)[ni].heap.replace(it->second, r);
deba@326
   975
                  (*_node_data)[ni].heap.decrease(r, rw);
deba@326
   976
                  it->second = r;
deba@326
   977
                }
deba@326
   978
              } else {
deba@326
   979
                (*_node_data)[ni].heap.push(r, rw);
deba@326
   980
                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
deba@326
   981
              }
deba@326
   982
deba@326
   983
              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
deba@326
   984
                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
deba@326
   985
deba@326
   986
                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
deba@326
   987
                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
deba@326
   988
                               (*_blossom_data)[blossom].offset);
deba@326
   989
                } else if ((*_delta2)[blossom] >
deba@326
   990
                           _blossom_set->classPrio(blossom) -
deba@326
   991
                           (*_blossom_data)[blossom].offset){
deba@326
   992
                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
deba@326
   993
                                   (*_blossom_data)[blossom].offset);
deba@326
   994
                }
deba@326
   995
              }
deba@326
   996
            }
deba@326
   997
deba@326
   998
          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
deba@326
   999
            if (_delta3->state(e) == _delta3->IN_HEAP) {
deba@326
  1000
              _delta3->erase(e);
deba@326
  1001
            }
deba@326
  1002
          } else {
deba@326
  1003
deba@326
  1004
            typename std::map<int, Arc>::iterator it =
deba@326
  1005
              (*_node_data)[vi].heap_index.find(tree);
deba@326
  1006
deba@326
  1007
            if (it != (*_node_data)[vi].heap_index.end()) {
deba@326
  1008
              (*_node_data)[vi].heap.erase(it->second);
deba@326
  1009
              (*_node_data)[vi].heap_index.erase(it);
deba@326
  1010
              if ((*_node_data)[vi].heap.empty()) {
deba@326
  1011
                _blossom_set->increase(v, std::numeric_limits<Value>::max());
deba@326
  1012
              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
deba@326
  1013
                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
deba@326
  1014
              }
deba@326
  1015
deba@326
  1016
              if ((*_blossom_data)[vb].status == MATCHED) {
deba@326
  1017
                if (_blossom_set->classPrio(vb) ==
deba@326
  1018
                    std::numeric_limits<Value>::max()) {
deba@326
  1019
                  _delta2->erase(vb);
deba@326
  1020
                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
deba@326
  1021
                           (*_blossom_data)[vb].offset) {
deba@326
  1022
                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
deba@326
  1023
                                   (*_blossom_data)[vb].offset);
deba@326
  1024
                }
deba@326
  1025
              }
deba@326
  1026
            }
deba@326
  1027
          }
deba@326
  1028
        }
deba@326
  1029
      }
deba@326
  1030
    }
deba@326
  1031
deba@326
  1032
    void oddToMatched(int blossom) {
deba@326
  1033
      (*_blossom_data)[blossom].offset -= _delta_sum;
deba@326
  1034
deba@326
  1035
      if (_blossom_set->classPrio(blossom) !=
deba@326
  1036
          std::numeric_limits<Value>::max()) {
deba@326
  1037
        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
deba@326
  1038
                       (*_blossom_data)[blossom].offset);
deba@326
  1039
      }
deba@326
  1040
deba@326
  1041
      if (!_blossom_set->trivial(blossom)) {
deba@326
  1042
        _delta4->erase(blossom);
deba@326
  1043
      }
deba@326
  1044
    }
deba@326
  1045
deba@326
  1046
    void oddToEven(int blossom, int tree) {
deba@326
  1047
      if (!_blossom_set->trivial(blossom)) {
deba@326
  1048
        _delta4->erase(blossom);
deba@326
  1049
        (*_blossom_data)[blossom].pot -=
deba@326
  1050
          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
deba@326
  1051
      }
deba@326
  1052
deba@326
  1053
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
deba@326
  1054
           n != INVALID; ++n) {
deba@326
  1055
        int ni = (*_node_index)[n];
deba@326
  1056
deba@326
  1057
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
deba@326
  1058
deba@326
  1059
        (*_node_data)[ni].heap.clear();
deba@326
  1060
        (*_node_data)[ni].heap_index.clear();
deba@326
  1061
        (*_node_data)[ni].pot +=
deba@326
  1062
          2 * _delta_sum - (*_blossom_data)[blossom].offset;
deba@326
  1063
deba@326
  1064
        _delta1->push(n, (*_node_data)[ni].pot);
deba@326
  1065
deba@326
  1066
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
deba@326
  1067
          Node v = _graph.source(e);
deba@326
  1068
          int vb = _blossom_set->find(v);
deba@326
  1069
          int vi = (*_node_index)[v];
deba@326
  1070
deba@326
  1071
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
deba@326
  1072
            dualScale * _weight[e];
deba@326
  1073
deba@326
  1074
          if ((*_blossom_data)[vb].status == EVEN) {
deba@326
  1075
            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
deba@326
  1076
              _delta3->push(e, rw / 2);
deba@326
  1077
            }
deba@326
  1078
          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
deba@326
  1079
            if (_delta3->state(e) != _delta3->IN_HEAP) {
deba@326
  1080
              _delta3->push(e, rw);
deba@326
  1081
            }
deba@326
  1082
          } else {
deba@326
  1083
deba@326
  1084
            typename std::map<int, Arc>::iterator it =
deba@326
  1085
              (*_node_data)[vi].heap_index.find(tree);
deba@326
  1086
deba@326
  1087
            if (it != (*_node_data)[vi].heap_index.end()) {
deba@326
  1088
              if ((*_node_data)[vi].heap[it->second] > rw) {
deba@326
  1089
                (*_node_data)[vi].heap.replace(it->second, e);
deba@326
  1090
                (*_node_data)[vi].heap.decrease(e, rw);
deba@326
  1091
                it->second = e;
deba@326
  1092
              }
deba@326
  1093
            } else {
deba@326
  1094
              (*_node_data)[vi].heap.push(e, rw);
deba@326
  1095
              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
deba@326
  1096
            }
deba@326
  1097
deba@326
  1098
            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
deba@326
  1099
              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
deba@326
  1100
deba@326
  1101
              if ((*_blossom_data)[vb].status == MATCHED) {
deba@326
  1102
                if (_delta2->state(vb) != _delta2->IN_HEAP) {
deba@326
  1103
                  _delta2->push(vb, _blossom_set->classPrio(vb) -
deba@326
  1104
                               (*_blossom_data)[vb].offset);
deba@326
  1105
                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
deba@326
  1106
                           (*_blossom_data)[vb].offset) {
deba@326
  1107
                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
deba@326
  1108
                                   (*_blossom_data)[vb].offset);
deba@326
  1109
                }
deba@326
  1110
              }
deba@326
  1111
            }
deba@326
  1112
          }
deba@326
  1113
        }
deba@326
  1114
      }
deba@326
  1115
      (*_blossom_data)[blossom].offset = 0;
deba@326
  1116
    }
deba@326
  1117
deba@326
  1118
deba@326
  1119
    void matchedToUnmatched(int blossom) {
deba@326
  1120
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
deba@326
  1121
        _delta2->erase(blossom);
deba@326
  1122
      }
deba@326
  1123
deba@326
  1124
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
deba@326
  1125
           n != INVALID; ++n) {
deba@326
  1126
        int ni = (*_node_index)[n];
deba@326
  1127
deba@326
  1128
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
deba@326
  1129
deba@326
  1130
        (*_node_data)[ni].heap.clear();
deba@326
  1131
        (*_node_data)[ni].heap_index.clear();
deba@326
  1132
deba@326
  1133
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
deba@326
  1134
          Node v = _graph.target(e);
deba@326
  1135
          int vb = _blossom_set->find(v);
deba@326
  1136
          int vi = (*_node_index)[v];
deba@326
  1137
deba@326
  1138
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
deba@326
  1139
            dualScale * _weight[e];
deba@326
  1140
deba@326
  1141
          if ((*_blossom_data)[vb].status == EVEN) {
deba@326
  1142
            if (_delta3->state(e) != _delta3->IN_HEAP) {
deba@326
  1143
              _delta3->push(e, rw);
deba@326
  1144
            }
deba@326
  1145
          }
deba@326
  1146
        }
deba@326
  1147
      }
deba@326
  1148
    }
deba@326
  1149
deba@326
  1150
    void unmatchedToMatched(int blossom) {
deba@326
  1151
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
deba@326
  1152
           n != INVALID; ++n) {
deba@326
  1153
        int ni = (*_node_index)[n];
deba@326
  1154
deba@326
  1155
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
deba@326
  1156
          Node v = _graph.source(e);
deba@326
  1157
          int vb = _blossom_set->find(v);
deba@326
  1158
          int vi = (*_node_index)[v];
deba@326
  1159
deba@326
  1160
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
deba@326
  1161
            dualScale * _weight[e];
deba@326
  1162
deba@326
  1163
          if (vb == blossom) {
deba@326
  1164
            if (_delta3->state(e) == _delta3->IN_HEAP) {
deba@326
  1165
              _delta3->erase(e);
deba@326
  1166
            }
deba@326
  1167
          } else if ((*_blossom_data)[vb].status == EVEN) {
deba@326
  1168
deba@326
  1169
            if (_delta3->state(e) == _delta3->IN_HEAP) {
deba@326
  1170
              _delta3->erase(e);
deba@326
  1171
            }
deba@326
  1172
deba@326
  1173
            int vt = _tree_set->find(vb);
deba@326
  1174
deba@326
  1175
            Arc r = _graph.oppositeArc(e);
deba@326
  1176
deba@326
  1177
            typename std::map<int, Arc>::iterator it =
deba@326
  1178
              (*_node_data)[ni].heap_index.find(vt);
deba@326
  1179
deba@326
  1180
            if (it != (*_node_data)[ni].heap_index.end()) {
deba@326
  1181
              if ((*_node_data)[ni].heap[it->second] > rw) {
deba@326
  1182
                (*_node_data)[ni].heap.replace(it->second, r);
deba@326
  1183
                (*_node_data)[ni].heap.decrease(r, rw);
deba@326
  1184
                it->second = r;
deba@326
  1185
              }
deba@326
  1186
            } else {
deba@326
  1187
              (*_node_data)[ni].heap.push(r, rw);
deba@326
  1188
              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
deba@326
  1189
            }
deba@326
  1190
deba@326
  1191
            if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
deba@326
  1192
              _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
deba@326
  1193
deba@326
  1194
              if (_delta2->state(blossom) != _delta2->IN_HEAP) {
deba@326
  1195
                _delta2->push(blossom, _blossom_set->classPrio(blossom) -
deba@326
  1196
                             (*_blossom_data)[blossom].offset);
deba@326
  1197
              } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
deba@326
  1198
                         (*_blossom_data)[blossom].offset){
deba@326
  1199
                _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
deba@326
  1200
                                 (*_blossom_data)[blossom].offset);
deba@326
  1201
              }
deba@326
  1202
            }
deba@326
  1203
deba@326
  1204
          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
deba@326
  1205
            if (_delta3->state(e) == _delta3->IN_HEAP) {
deba@326
  1206
              _delta3->erase(e);
deba@326
  1207
            }
deba@326
  1208
          }
deba@326
  1209
        }
deba@326
  1210
      }
deba@326
  1211
    }
deba@326
  1212
deba@326
  1213
    void alternatePath(int even, int tree) {
deba@326
  1214
      int odd;
deba@326
  1215
deba@326
  1216
      evenToMatched(even, tree);
deba@326
  1217
      (*_blossom_data)[even].status = MATCHED;
deba@326
  1218
deba@326
  1219
      while ((*_blossom_data)[even].pred != INVALID) {
deba@326
  1220
        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
deba@326
  1221
        (*_blossom_data)[odd].status = MATCHED;
deba@326
  1222
        oddToMatched(odd);
deba@326
  1223
        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
deba@326
  1224
deba@326
  1225
        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
deba@326
  1226
        (*_blossom_data)[even].status = MATCHED;
deba@326
  1227
        evenToMatched(even, tree);
deba@326
  1228
        (*_blossom_data)[even].next =
deba@326
  1229
          _graph.oppositeArc((*_blossom_data)[odd].pred);
deba@326
  1230
      }
deba@326
  1231
deba@326
  1232
    }
deba@326
  1233
deba@326
  1234
    void destroyTree(int tree) {
deba@326
  1235
      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
deba@326
  1236
        if ((*_blossom_data)[b].status == EVEN) {
deba@326
  1237
          (*_blossom_data)[b].status = MATCHED;
deba@326
  1238
          evenToMatched(b, tree);
deba@326
  1239
        } else if ((*_blossom_data)[b].status == ODD) {
deba@326
  1240
          (*_blossom_data)[b].status = MATCHED;
deba@326
  1241
          oddToMatched(b);
deba@326
  1242
        }
deba@326
  1243
      }
deba@326
  1244
      _tree_set->eraseClass(tree);
deba@326
  1245
    }
deba@326
  1246
deba@326
  1247
deba@326
  1248
    void unmatchNode(const Node& node) {
deba@326
  1249
      int blossom = _blossom_set->find(node);
deba@326
  1250
      int tree = _tree_set->find(blossom);
deba@326
  1251
deba@326
  1252
      alternatePath(blossom, tree);
deba@326
  1253
      destroyTree(tree);
deba@326
  1254
deba@326
  1255
      (*_blossom_data)[blossom].status = UNMATCHED;
deba@326
  1256
      (*_blossom_data)[blossom].base = node;
deba@326
  1257
      matchedToUnmatched(blossom);
deba@326
  1258
    }
deba@326
  1259
deba@326
  1260
deba@327
  1261
    void augmentOnEdge(const Edge& edge) {
deba@327
  1262
deba@327
  1263
      int left = _blossom_set->find(_graph.u(edge));
deba@327
  1264
      int right = _blossom_set->find(_graph.v(edge));
deba@326
  1265
deba@326
  1266
      if ((*_blossom_data)[left].status == EVEN) {
deba@326
  1267
        int left_tree = _tree_set->find(left);
deba@326
  1268
        alternatePath(left, left_tree);
deba@326
  1269
        destroyTree(left_tree);
deba@326
  1270
      } else {
deba@326
  1271
        (*_blossom_data)[left].status = MATCHED;
deba@326
  1272
        unmatchedToMatched(left);
deba@326
  1273
      }
deba@326
  1274
deba@326
  1275
      if ((*_blossom_data)[right].status == EVEN) {
deba@326
  1276
        int right_tree = _tree_set->find(right);
deba@326
  1277
        alternatePath(right, right_tree);
deba@326
  1278
        destroyTree(right_tree);
deba@326
  1279
      } else {
deba@326
  1280
        (*_blossom_data)[right].status = MATCHED;
deba@326
  1281
        unmatchedToMatched(right);
deba@326
  1282
      }
deba@326
  1283
deba@327
  1284
      (*_blossom_data)[left].next = _graph.direct(edge, true);
deba@327
  1285
      (*_blossom_data)[right].next = _graph.direct(edge, false);
deba@326
  1286
    }
deba@326
  1287
deba@326
  1288
    void extendOnArc(const Arc& arc) {
deba@326
  1289
      int base = _blossom_set->find(_graph.target(arc));
deba@326
  1290
      int tree = _tree_set->find(base);
deba@326
  1291
deba@326
  1292
      int odd = _blossom_set->find(_graph.source(arc));
deba@326
  1293
      _tree_set->insert(odd, tree);
deba@326
  1294
      (*_blossom_data)[odd].status = ODD;
deba@326
  1295
      matchedToOdd(odd);
deba@326
  1296
      (*_blossom_data)[odd].pred = arc;
deba@326
  1297
deba@326
  1298
      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
deba@326
  1299
      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
deba@326
  1300
      _tree_set->insert(even, tree);
deba@326
  1301
      (*_blossom_data)[even].status = EVEN;
deba@326
  1302
      matchedToEven(even, tree);
deba@326
  1303
    }
deba@326
  1304
deba@327
  1305
    void shrinkOnEdge(const Edge& edge, int tree) {
deba@326
  1306
      int nca = -1;
deba@326
  1307
      std::vector<int> left_path, right_path;
deba@326
  1308
deba@326
  1309
      {
deba@326
  1310
        std::set<int> left_set, right_set;
deba@326
  1311
        int left = _blossom_set->find(_graph.u(edge));
deba@326
  1312
        left_path.push_back(left);
deba@326
  1313
        left_set.insert(left);
deba@326
  1314
deba@326
  1315
        int right = _blossom_set->find(_graph.v(edge));
deba@326
  1316
        right_path.push_back(right);
deba@326
  1317
        right_set.insert(right);
deba@326
  1318
deba@326
  1319
        while (true) {
deba@326
  1320
deba@326
  1321
          if ((*_blossom_data)[left].pred == INVALID) break;
deba@326
  1322
deba@326
  1323
          left =
deba@326
  1324
            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
deba@326
  1325
          left_path.push_back(left);
deba@326
  1326
          left =
deba@326
  1327
            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
deba@326
  1328
          left_path.push_back(left);
deba@326
  1329
deba@326
  1330
          left_set.insert(left);
deba@326
  1331
deba@326
  1332
          if (right_set.find(left) != right_set.end()) {
deba@326
  1333
            nca = left;
deba@326
  1334
            break;
deba@326
  1335
          }
deba@326
  1336
deba@326
  1337
          if ((*_blossom_data)[right].pred == INVALID) break;
deba@326
  1338
deba@326
  1339
          right =
deba@326
  1340
            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
deba@326
  1341
          right_path.push_back(right);
deba@326
  1342
          right =
deba@326
  1343
            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
deba@326
  1344
          right_path.push_back(right);
deba@326
  1345
deba@326
  1346
          right_set.insert(right);
deba@326
  1347
deba@326
  1348
          if (left_set.find(right) != left_set.end()) {
deba@326
  1349
            nca = right;
deba@326
  1350
            break;
deba@326
  1351
          }
deba@326
  1352
deba@326
  1353
        }
deba@326
  1354
deba@326
  1355
        if (nca == -1) {
deba@326
  1356
          if ((*_blossom_data)[left].pred == INVALID) {
deba@326
  1357
            nca = right;
deba@326
  1358
            while (left_set.find(nca) == left_set.end()) {
deba@326
  1359
              nca =
deba@326
  1360
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
deba@326
  1361
              right_path.push_back(nca);
deba@326
  1362
              nca =
deba@326
  1363
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
deba@326
  1364
              right_path.push_back(nca);
deba@326
  1365
            }
deba@326
  1366
          } else {
deba@326
  1367
            nca = left;
deba@326
  1368
            while (right_set.find(nca) == right_set.end()) {
deba@326
  1369
              nca =
deba@326
  1370
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
deba@326
  1371
              left_path.push_back(nca);
deba@326
  1372
              nca =
deba@326
  1373
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
deba@326
  1374
              left_path.push_back(nca);
deba@326
  1375
            }
deba@326
  1376
          }
deba@326
  1377
        }
deba@326
  1378
      }
deba@326
  1379
deba@326
  1380
      std::vector<int> subblossoms;
deba@326
  1381
      Arc prev;
deba@326
  1382
deba@326
  1383
      prev = _graph.direct(edge, true);
deba@326
  1384
      for (int i = 0; left_path[i] != nca; i += 2) {
deba@326
  1385
        subblossoms.push_back(left_path[i]);
deba@326
  1386
        (*_blossom_data)[left_path[i]].next = prev;
deba@326
  1387
        _tree_set->erase(left_path[i]);
deba@326
  1388
deba@326
  1389
        subblossoms.push_back(left_path[i + 1]);
deba@326
  1390
        (*_blossom_data)[left_path[i + 1]].status = EVEN;
deba@326
  1391
        oddToEven(left_path[i + 1], tree);
deba@326
  1392
        _tree_set->erase(left_path[i + 1]);
deba@326
  1393
        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
deba@326
  1394
      }
deba@326
  1395
deba@326
  1396
      int k = 0;
deba@326
  1397
      while (right_path[k] != nca) ++k;
deba@326
  1398
deba@326
  1399
      subblossoms.push_back(nca);
deba@326
  1400
      (*_blossom_data)[nca].next = prev;
deba@326
  1401
deba@326
  1402
      for (int i = k - 2; i >= 0; i -= 2) {
deba@326
  1403
        subblossoms.push_back(right_path[i + 1]);
deba@326
  1404
        (*_blossom_data)[right_path[i + 1]].status = EVEN;
deba@326
  1405
        oddToEven(right_path[i + 1], tree);
deba@326
  1406
        _tree_set->erase(right_path[i + 1]);
deba@326
  1407
deba@326
  1408
        (*_blossom_data)[right_path[i + 1]].next =
deba@326
  1409
          (*_blossom_data)[right_path[i + 1]].pred;
deba@326
  1410
deba@326
  1411
        subblossoms.push_back(right_path[i]);
deba@326
  1412
        _tree_set->erase(right_path[i]);
deba@326
  1413
      }
deba@326
  1414
deba@326
  1415
      int surface =
deba@326
  1416
        _blossom_set->join(subblossoms.begin(), subblossoms.end());
deba@326
  1417
deba@326
  1418
      for (int i = 0; i < int(subblossoms.size()); ++i) {
deba@326
  1419
        if (!_blossom_set->trivial(subblossoms[i])) {
deba@326
  1420
          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
deba@326
  1421
        }
deba@326
  1422
        (*_blossom_data)[subblossoms[i]].status = MATCHED;
deba@326
  1423
      }
deba@326
  1424
deba@326
  1425
      (*_blossom_data)[surface].pot = -2 * _delta_sum;
deba@326
  1426
      (*_blossom_data)[surface].offset = 0;
deba@326
  1427
      (*_blossom_data)[surface].status = EVEN;
deba@326
  1428
      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
deba@326
  1429
      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
deba@326
  1430
deba@326
  1431
      _tree_set->insert(surface, tree);
deba@326
  1432
      _tree_set->erase(nca);
deba@326
  1433
    }
deba@326
  1434
deba@326
  1435
    void splitBlossom(int blossom) {
deba@326
  1436
      Arc next = (*_blossom_data)[blossom].next;
deba@326
  1437
      Arc pred = (*_blossom_data)[blossom].pred;
deba@326
  1438
deba@326
  1439
      int tree = _tree_set->find(blossom);
deba@326
  1440
deba@326
  1441
      (*_blossom_data)[blossom].status = MATCHED;
deba@326
  1442
      oddToMatched(blossom);
deba@326
  1443
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
deba@326
  1444
        _delta2->erase(blossom);
deba@326
  1445
      }
deba@326
  1446
deba@326
  1447
      std::vector<int> subblossoms;
deba@326
  1448
      _blossom_set->split(blossom, std::back_inserter(subblossoms));
deba@326
  1449
deba@326
  1450
      Value offset = (*_blossom_data)[blossom].offset;
deba@326
  1451
      int b = _blossom_set->find(_graph.source(pred));
deba@326
  1452
      int d = _blossom_set->find(_graph.source(next));
deba@326
  1453
deba@326
  1454
      int ib = -1, id = -1;
deba@326
  1455
      for (int i = 0; i < int(subblossoms.size()); ++i) {
deba@326
  1456
        if (subblossoms[i] == b) ib = i;
deba@326
  1457
        if (subblossoms[i] == d) id = i;
deba@326
  1458
deba@326
  1459
        (*_blossom_data)[subblossoms[i]].offset = offset;
deba@326
  1460
        if (!_blossom_set->trivial(subblossoms[i])) {
deba@326
  1461
          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
deba@326
  1462
        }
deba@326
  1463
        if (_blossom_set->classPrio(subblossoms[i]) !=
deba@326
  1464
            std::numeric_limits<Value>::max()) {
deba@326
  1465
          _delta2->push(subblossoms[i],
deba@326
  1466
                        _blossom_set->classPrio(subblossoms[i]) -
deba@326
  1467
                        (*_blossom_data)[subblossoms[i]].offset);
deba@326
  1468
        }
deba@326
  1469
      }
deba@326
  1470
deba@326
  1471
      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
deba@326
  1472
        for (int i = (id + 1) % subblossoms.size();
deba@326
  1473
             i != ib; i = (i + 2) % subblossoms.size()) {
deba@326
  1474
          int sb = subblossoms[i];
deba@326
  1475
          int tb = subblossoms[(i + 1) % subblossoms.size()];
deba@326
  1476
          (*_blossom_data)[sb].next =
deba@326
  1477
            _graph.oppositeArc((*_blossom_data)[tb].next);
deba@326
  1478
        }
deba@326
  1479
deba@326
  1480
        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
deba@326
  1481
          int sb = subblossoms[i];
deba@326
  1482
          int tb = subblossoms[(i + 1) % subblossoms.size()];
deba@326
  1483
          int ub = subblossoms[(i + 2) % subblossoms.size()];
deba@326
  1484
deba@326
  1485
          (*_blossom_data)[sb].status = ODD;
deba@326
  1486
          matchedToOdd(sb);
deba@326
  1487
          _tree_set->insert(sb, tree);
deba@326
  1488
          (*_blossom_data)[sb].pred = pred;
deba@326
  1489
          (*_blossom_data)[sb].next =
deba@326
  1490
                           _graph.oppositeArc((*_blossom_data)[tb].next);
deba@326
  1491
deba@326
  1492
          pred = (*_blossom_data)[ub].next;
deba@326
  1493
deba@326
  1494
          (*_blossom_data)[tb].status = EVEN;
deba@326
  1495
          matchedToEven(tb, tree);
deba@326
  1496
          _tree_set->insert(tb, tree);
deba@326
  1497
          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
deba@326
  1498
        }
deba@326
  1499
deba@326
  1500
        (*_blossom_data)[subblossoms[id]].status = ODD;
deba@326
  1501
        matchedToOdd(subblossoms[id]);
deba@326
  1502
        _tree_set->insert(subblossoms[id], tree);
deba@326
  1503
        (*_blossom_data)[subblossoms[id]].next = next;
deba@326
  1504
        (*_blossom_data)[subblossoms[id]].pred = pred;
deba@326
  1505
deba@326
  1506
      } else {
deba@326
  1507
deba@326
  1508
        for (int i = (ib + 1) % subblossoms.size();
deba@326
  1509
             i != id; i = (i + 2) % subblossoms.size()) {
deba@326
  1510
          int sb = subblossoms[i];
deba@326
  1511
          int tb = subblossoms[(i + 1) % subblossoms.size()];
deba@326
  1512
          (*_blossom_data)[sb].next =
deba@326
  1513
            _graph.oppositeArc((*_blossom_data)[tb].next);
deba@326
  1514
        }
deba@326
  1515
deba@326
  1516
        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
deba@326
  1517
          int sb = subblossoms[i];
deba@326
  1518
          int tb = subblossoms[(i + 1) % subblossoms.size()];
deba@326
  1519
          int ub = subblossoms[(i + 2) % subblossoms.size()];
deba@326
  1520
deba@326
  1521
          (*_blossom_data)[sb].status = ODD;
deba@326
  1522
          matchedToOdd(sb);
deba@326
  1523
          _tree_set->insert(sb, tree);
deba@326
  1524
          (*_blossom_data)[sb].next = next;
deba@326
  1525
          (*_blossom_data)[sb].pred =
deba@326
  1526
            _graph.oppositeArc((*_blossom_data)[tb].next);
deba@326
  1527
deba@326
  1528
          (*_blossom_data)[tb].status = EVEN;
deba@326
  1529
          matchedToEven(tb, tree);
deba@326
  1530
          _tree_set->insert(tb, tree);
deba@326
  1531
          (*_blossom_data)[tb].pred =
deba@326
  1532
            (*_blossom_data)[tb].next =
deba@326
  1533
            _graph.oppositeArc((*_blossom_data)[ub].next);
deba@326
  1534
          next = (*_blossom_data)[ub].next;
deba@326
  1535
        }
deba@326
  1536
deba@326
  1537
        (*_blossom_data)[subblossoms[ib]].status = ODD;
deba@326
  1538
        matchedToOdd(subblossoms[ib]);
deba@326
  1539
        _tree_set->insert(subblossoms[ib], tree);
deba@326
  1540
        (*_blossom_data)[subblossoms[ib]].next = next;
deba@326
  1541
        (*_blossom_data)[subblossoms[ib]].pred = pred;
deba@326
  1542
      }
deba@326
  1543
      _tree_set->erase(blossom);
deba@326
  1544
    }
deba@326
  1545
deba@326
  1546
    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
deba@326
  1547
      if (_blossom_set->trivial(blossom)) {
deba@326
  1548
        int bi = (*_node_index)[base];
deba@326
  1549
        Value pot = (*_node_data)[bi].pot;
deba@326
  1550
kpeter@581
  1551
        (*_matching)[base] = matching;
deba@326
  1552
        _blossom_node_list.push_back(base);
kpeter@581
  1553
        (*_node_potential)[base] = pot;
deba@326
  1554
      } else {
deba@326
  1555
deba@326
  1556
        Value pot = (*_blossom_data)[blossom].pot;
deba@326
  1557
        int bn = _blossom_node_list.size();
deba@326
  1558
deba@326
  1559
        std::vector<int> subblossoms;
deba@326
  1560
        _blossom_set->split(blossom, std::back_inserter(subblossoms));
deba@326
  1561
        int b = _blossom_set->find(base);
deba@326
  1562
        int ib = -1;
deba@326
  1563
        for (int i = 0; i < int(subblossoms.size()); ++i) {
deba@326
  1564
          if (subblossoms[i] == b) { ib = i; break; }
deba@326
  1565
        }
deba@326
  1566
deba@326
  1567
        for (int i = 1; i < int(subblossoms.size()); i += 2) {
deba@326
  1568
          int sb = subblossoms[(ib + i) % subblossoms.size()];
deba@326
  1569
          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
deba@326
  1570
deba@326
  1571
          Arc m = (*_blossom_data)[tb].next;
deba@326
  1572
          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
deba@326
  1573
          extractBlossom(tb, _graph.source(m), m);
deba@326
  1574
        }
deba@326
  1575
        extractBlossom(subblossoms[ib], base, matching);
deba@326
  1576
deba@326
  1577
        int en = _blossom_node_list.size();
deba@326
  1578
deba@326
  1579
        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
deba@326
  1580
      }
deba@326
  1581
    }
deba@326
  1582
deba@326
  1583
    void extractMatching() {
deba@326
  1584
      std::vector<int> blossoms;
deba@326
  1585
      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
deba@326
  1586
        blossoms.push_back(c);
deba@326
  1587
      }
deba@326
  1588
deba@326
  1589
      for (int i = 0; i < int(blossoms.size()); ++i) {
deba@326
  1590
        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
deba@326
  1591
deba@326
  1592
          Value offset = (*_blossom_data)[blossoms[i]].offset;
deba@326
  1593
          (*_blossom_data)[blossoms[i]].pot += 2 * offset;
deba@326
  1594
          for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
deba@326
  1595
               n != INVALID; ++n) {
deba@326
  1596
            (*_node_data)[(*_node_index)[n]].pot -= offset;
deba@326
  1597
          }
deba@326
  1598
deba@326
  1599
          Arc matching = (*_blossom_data)[blossoms[i]].next;
deba@326
  1600
          Node base = _graph.source(matching);
deba@326
  1601
          extractBlossom(blossoms[i], base, matching);
deba@326
  1602
        } else {
deba@326
  1603
          Node base = (*_blossom_data)[blossoms[i]].base;
deba@326
  1604
          extractBlossom(blossoms[i], base, INVALID);
deba@326
  1605
        }
deba@326
  1606
      }
deba@326
  1607
    }
deba@326
  1608
deba@326
  1609
  public:
deba@326
  1610
deba@326
  1611
    /// \brief Constructor
deba@326
  1612
    ///
deba@326
  1613
    /// Constructor.
deba@326
  1614
    MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
deba@326
  1615
      : _graph(graph), _weight(weight), _matching(0),
deba@326
  1616
        _node_potential(0), _blossom_potential(), _blossom_node_list(),
deba@326
  1617
        _node_num(0), _blossom_num(0),
deba@326
  1618
deba@326
  1619
        _blossom_index(0), _blossom_set(0), _blossom_data(0),
deba@326
  1620
        _node_index(0), _node_heap_index(0), _node_data(0),
deba@326
  1621
        _tree_set_index(0), _tree_set(0),
deba@326
  1622
deba@326
  1623
        _delta1_index(0), _delta1(0),
deba@326
  1624
        _delta2_index(0), _delta2(0),
deba@326
  1625
        _delta3_index(0), _delta3(0),
deba@326
  1626
        _delta4_index(0), _delta4(0),
deba@326
  1627
deba@326
  1628
        _delta_sum() {}
deba@326
  1629
deba@326
  1630
    ~MaxWeightedMatching() {
deba@326
  1631
      destroyStructures();
deba@326
  1632
    }
deba@326
  1633
deba@326
  1634
    /// \name Execution control
alpar@330
  1635
    /// The simplest way to execute the algorithm is to use the
deba@326
  1636
    /// \c run() member function.
deba@326
  1637
deba@326
  1638
    ///@{
deba@326
  1639
deba@326
  1640
    /// \brief Initialize the algorithm
deba@326
  1641
    ///
deba@326
  1642
    /// Initialize the algorithm
deba@326
  1643
    void init() {
deba@326
  1644
      createStructures();
deba@326
  1645
deba@326
  1646
      for (ArcIt e(_graph); e != INVALID; ++e) {
kpeter@581
  1647
        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
deba@326
  1648
      }
deba@326
  1649
      for (NodeIt n(_graph); n != INVALID; ++n) {
kpeter@581
  1650
        (*_delta1_index)[n] = _delta1->PRE_HEAP;
deba@326
  1651
      }
deba@326
  1652
      for (EdgeIt e(_graph); e != INVALID; ++e) {
kpeter@581
  1653
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
deba@326
  1654
      }
deba@326
  1655
      for (int i = 0; i < _blossom_num; ++i) {
kpeter@581
  1656
        (*_delta2_index)[i] = _delta2->PRE_HEAP;
kpeter@581
  1657
        (*_delta4_index)[i] = _delta4->PRE_HEAP;
deba@326
  1658
      }
deba@326
  1659
deba@326
  1660
      int index = 0;
deba@326
  1661
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@326
  1662
        Value max = 0;
deba@326
  1663
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
deba@326
  1664
          if (_graph.target(e) == n) continue;
deba@326
  1665
          if ((dualScale * _weight[e]) / 2 > max) {
deba@326
  1666
            max = (dualScale * _weight[e]) / 2;
deba@326
  1667
          }
deba@326
  1668
        }
kpeter@581
  1669
        (*_node_index)[n] = index;
deba@326
  1670
        (*_node_data)[index].pot = max;
deba@326
  1671
        _delta1->push(n, max);
deba@326
  1672
        int blossom =
deba@326
  1673
          _blossom_set->insert(n, std::numeric_limits<Value>::max());
deba@326
  1674
deba@326
  1675
        _tree_set->insert(blossom);
deba@326
  1676
deba@326
  1677
        (*_blossom_data)[blossom].status = EVEN;
deba@326
  1678
        (*_blossom_data)[blossom].pred = INVALID;
deba@326
  1679
        (*_blossom_data)[blossom].next = INVALID;
deba@326
  1680
        (*_blossom_data)[blossom].pot = 0;
deba@326
  1681
        (*_blossom_data)[blossom].offset = 0;
deba@326
  1682
        ++index;
deba@326
  1683
      }
deba@326
  1684
      for (EdgeIt e(_graph); e != INVALID; ++e) {
deba@326
  1685
        int si = (*_node_index)[_graph.u(e)];
deba@326
  1686
        int ti = (*_node_index)[_graph.v(e)];
deba@326
  1687
        if (_graph.u(e) != _graph.v(e)) {
deba@326
  1688
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
deba@326
  1689
                            dualScale * _weight[e]) / 2);
deba@326
  1690
        }
deba@326
  1691
      }
deba@326
  1692
    }
deba@326
  1693
deba@326
  1694
    /// \brief Starts the algorithm
deba@326
  1695
    ///
deba@326
  1696
    /// Starts the algorithm
deba@326
  1697
    void start() {
deba@326
  1698
      enum OpType {
deba@326
  1699
        D1, D2, D3, D4
deba@326
  1700
      };
deba@326
  1701
deba@326
  1702
      int unmatched = _node_num;
deba@326
  1703
      while (unmatched > 0) {
deba@326
  1704
        Value d1 = !_delta1->empty() ?
deba@326
  1705
          _delta1->prio() : std::numeric_limits<Value>::max();
deba@326
  1706
deba@326
  1707
        Value d2 = !_delta2->empty() ?
deba@326
  1708
          _delta2->prio() : std::numeric_limits<Value>::max();
deba@326
  1709
deba@326
  1710
        Value d3 = !_delta3->empty() ?
deba@326
  1711
          _delta3->prio() : std::numeric_limits<Value>::max();
deba@326
  1712
deba@326
  1713
        Value d4 = !_delta4->empty() ?
deba@326
  1714
          _delta4->prio() : std::numeric_limits<Value>::max();
deba@326
  1715
deba@326
  1716
        _delta_sum = d1; OpType ot = D1;
deba@326
  1717
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
deba@326
  1718
        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
deba@326
  1719
        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
deba@326
  1720
deba@326
  1721
deba@326
  1722
        switch (ot) {
deba@326
  1723
        case D1:
deba@326
  1724
          {
deba@326
  1725
            Node n = _delta1->top();
deba@326
  1726
            unmatchNode(n);
deba@326
  1727
            --unmatched;
deba@326
  1728
          }
deba@326
  1729
          break;
deba@326
  1730
        case D2:
deba@326
  1731
          {
deba@326
  1732
            int blossom = _delta2->top();
deba@326
  1733
            Node n = _blossom_set->classTop(blossom);
deba@326
  1734
            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
deba@326
  1735
            extendOnArc(e);
deba@326
  1736
          }
deba@326
  1737
          break;
deba@326
  1738
        case D3:
deba@326
  1739
          {
deba@326
  1740
            Edge e = _delta3->top();
deba@326
  1741
deba@326
  1742
            int left_blossom = _blossom_set->find(_graph.u(e));
deba@326
  1743
            int right_blossom = _blossom_set->find(_graph.v(e));
deba@326
  1744
deba@326
  1745
            if (left_blossom == right_blossom) {
deba@326
  1746
              _delta3->pop();
deba@326
  1747
            } else {
deba@326
  1748
              int left_tree;
deba@326
  1749
              if ((*_blossom_data)[left_blossom].status == EVEN) {
deba@326
  1750
                left_tree = _tree_set->find(left_blossom);
deba@326
  1751
              } else {
deba@326
  1752
                left_tree = -1;
deba@326
  1753
                ++unmatched;
deba@326
  1754
              }
deba@326
  1755
              int right_tree;
deba@326
  1756
              if ((*_blossom_data)[right_blossom].status == EVEN) {
deba@326
  1757
                right_tree = _tree_set->find(right_blossom);
deba@326
  1758
              } else {
deba@326
  1759
                right_tree = -1;
deba@326
  1760
                ++unmatched;
deba@326
  1761
              }
deba@326
  1762
deba@326
  1763
              if (left_tree == right_tree) {
deba@327
  1764
                shrinkOnEdge(e, left_tree);
deba@326
  1765
              } else {
deba@327
  1766
                augmentOnEdge(e);
deba@326
  1767
                unmatched -= 2;
deba@326
  1768
              }
deba@326
  1769
            }
deba@326
  1770
          } break;
deba@326
  1771
        case D4:
deba@326
  1772
          splitBlossom(_delta4->top());
deba@326
  1773
          break;
deba@326
  1774
        }
deba@326
  1775
      }
deba@326
  1776
      extractMatching();
deba@326
  1777
    }
deba@326
  1778
deba@326
  1779
    /// \brief Runs %MaxWeightedMatching algorithm.
deba@326
  1780
    ///
deba@326
  1781
    /// This method runs the %MaxWeightedMatching algorithm.
deba@326
  1782
    ///
deba@326
  1783
    /// \note mwm.run() is just a shortcut of the following code.
deba@326
  1784
    /// \code
deba@326
  1785
    ///   mwm.init();
deba@326
  1786
    ///   mwm.start();
deba@326
  1787
    /// \endcode
deba@326
  1788
    void run() {
deba@326
  1789
      init();
deba@326
  1790
      start();
deba@326
  1791
    }
deba@326
  1792
deba@326
  1793
    /// @}
deba@326
  1794
deba@326
  1795
    /// \name Primal solution
alpar@330
  1796
    /// Functions to get the primal solution, ie. the matching.
deba@326
  1797
deba@326
  1798
    /// @{
deba@326
  1799
alpar@330
  1800
    /// \brief Returns the weight of the matching.
deba@326
  1801
    ///
alpar@330
  1802
    /// Returns the weight of the matching.
deba@326
  1803
    Value matchingValue() const {
deba@326
  1804
      Value sum = 0;
deba@326
  1805
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@326
  1806
        if ((*_matching)[n] != INVALID) {
deba@326
  1807
          sum += _weight[(*_matching)[n]];
deba@326
  1808
        }
deba@326
  1809
      }
deba@326
  1810
      return sum /= 2;
deba@326
  1811
    }
deba@326
  1812
deba@327
  1813
    /// \brief Returns the cardinality of the matching.
deba@326
  1814
    ///
deba@327
  1815
    /// Returns the cardinality of the matching.
deba@327
  1816
    int matchingSize() const {
deba@327
  1817
      int num = 0;
deba@327
  1818
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@327
  1819
        if ((*_matching)[n] != INVALID) {
deba@327
  1820
          ++num;
deba@327
  1821
        }
deba@327
  1822
      }
deba@327
  1823
      return num /= 2;
deba@327
  1824
    }
deba@327
  1825
deba@327
  1826
    /// \brief Returns true when the edge is in the matching.
deba@327
  1827
    ///
deba@327
  1828
    /// Returns true when the edge is in the matching.
deba@327
  1829
    bool matching(const Edge& edge) const {
deba@327
  1830
      return edge == (*_matching)[_graph.u(edge)];
deba@326
  1831
    }
deba@326
  1832
deba@326
  1833
    /// \brief Returns the incident matching arc.
deba@326
  1834
    ///
deba@326
  1835
    /// Returns the incident matching arc from given node. If the
deba@326
  1836
    /// node is not matched then it gives back \c INVALID.
deba@326
  1837
    Arc matching(const Node& node) const {
deba@326
  1838
      return (*_matching)[node];
deba@326
  1839
    }
deba@326
  1840
deba@326
  1841
    /// \brief Returns the mate of the node.
deba@326
  1842
    ///
deba@326
  1843
    /// Returns the adjancent node in a mathcing arc. If the node is
deba@326
  1844
    /// not matched then it gives back \c INVALID.
deba@326
  1845
    Node mate(const Node& node) const {
deba@326
  1846
      return (*_matching)[node] != INVALID ?
deba@326
  1847
        _graph.target((*_matching)[node]) : INVALID;
deba@326
  1848
    }
deba@326
  1849
deba@326
  1850
    /// @}
deba@326
  1851
deba@326
  1852
    /// \name Dual solution
alpar@330
  1853
    /// Functions to get the dual solution.
deba@326
  1854
deba@326
  1855
    /// @{
deba@326
  1856
deba@326
  1857
    /// \brief Returns the value of the dual solution.
deba@326
  1858
    ///
deba@326
  1859
    /// Returns the value of the dual solution. It should be equal to
deba@326
  1860
    /// the primal value scaled by \ref dualScale "dual scale".
deba@326
  1861
    Value dualValue() const {
deba@326
  1862
      Value sum = 0;
deba@326
  1863
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@326
  1864
        sum += nodeValue(n);
deba@326
  1865
      }
deba@326
  1866
      for (int i = 0; i < blossomNum(); ++i) {
deba@326
  1867
        sum += blossomValue(i) * (blossomSize(i) / 2);
deba@326
  1868
      }
deba@326
  1869
      return sum;
deba@326
  1870
    }
deba@326
  1871
deba@326
  1872
    /// \brief Returns the value of the node.
deba@326
  1873
    ///
deba@326
  1874
    /// Returns the the value of the node.
deba@326
  1875
    Value nodeValue(const Node& n) const {
deba@326
  1876
      return (*_node_potential)[n];
deba@326
  1877
    }
deba@326
  1878
deba@326
  1879
    /// \brief Returns the number of the blossoms in the basis.
deba@326
  1880
    ///
deba@326
  1881
    /// Returns the number of the blossoms in the basis.
deba@326
  1882
    /// \see BlossomIt
deba@326
  1883
    int blossomNum() const {
deba@326
  1884
      return _blossom_potential.size();
deba@326
  1885
    }
deba@326
  1886
deba@326
  1887
deba@326
  1888
    /// \brief Returns the number of the nodes in the blossom.
deba@326
  1889
    ///
deba@326
  1890
    /// Returns the number of the nodes in the blossom.
deba@326
  1891
    int blossomSize(int k) const {
deba@326
  1892
      return _blossom_potential[k].end - _blossom_potential[k].begin;
deba@326
  1893
    }
deba@326
  1894
deba@326
  1895
    /// \brief Returns the value of the blossom.
deba@326
  1896
    ///
deba@326
  1897
    /// Returns the the value of the blossom.
deba@326
  1898
    /// \see BlossomIt
deba@326
  1899
    Value blossomValue(int k) const {
deba@326
  1900
      return _blossom_potential[k].value;
deba@326
  1901
    }
deba@326
  1902
alpar@330
  1903
    /// \brief Iterator for obtaining the nodes of the blossom.
deba@326
  1904
    ///
alpar@330
  1905
    /// Iterator for obtaining the nodes of the blossom. This class
alpar@330
  1906
    /// provides a common lemon style iterator for listing a
deba@326
  1907
    /// subset of the nodes.
deba@326
  1908
    class BlossomIt {
deba@326
  1909
    public:
deba@326
  1910
deba@326
  1911
      /// \brief Constructor.
deba@326
  1912
      ///
alpar@330
  1913
      /// Constructor to get the nodes of the variable.
deba@326
  1914
      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
deba@326
  1915
        : _algorithm(&algorithm)
deba@326
  1916
      {
deba@326
  1917
        _index = _algorithm->_blossom_potential[variable].begin;
deba@326
  1918
        _last = _algorithm->_blossom_potential[variable].end;
deba@326
  1919
      }
deba@326
  1920
deba@326
  1921
      /// \brief Conversion to node.
deba@326
  1922
      ///
deba@326
  1923
      /// Conversion to node.
deba@326
  1924
      operator Node() const {
deba@327
  1925
        return _algorithm->_blossom_node_list[_index];
deba@326
  1926
      }
deba@326
  1927
deba@326
  1928
      /// \brief Increment operator.
deba@326
  1929
      ///
deba@326
  1930
      /// Increment operator.
deba@326
  1931
      BlossomIt& operator++() {
deba@326
  1932
        ++_index;
deba@326
  1933
        return *this;
deba@326
  1934
      }
deba@326
  1935
deba@327
  1936
      /// \brief Validity checking
deba@327
  1937
      ///
deba@327
  1938
      /// Checks whether the iterator is invalid.
deba@327
  1939
      bool operator==(Invalid) const { return _index == _last; }
deba@327
  1940
deba@327
  1941
      /// \brief Validity checking
deba@327
  1942
      ///
deba@327
  1943
      /// Checks whether the iterator is valid.
deba@327
  1944
      bool operator!=(Invalid) const { return _index != _last; }
deba@326
  1945
deba@326
  1946
    private:
deba@326
  1947
      const MaxWeightedMatching* _algorithm;
deba@326
  1948
      int _last;
deba@326
  1949
      int _index;
deba@326
  1950
    };
deba@326
  1951
deba@326
  1952
    /// @}
deba@326
  1953
deba@326
  1954
  };
deba@326
  1955
deba@326
  1956
  /// \ingroup matching
deba@326
  1957
  ///
deba@326
  1958
  /// \brief Weighted perfect matching in general graphs
deba@326
  1959
  ///
deba@326
  1960
  /// This class provides an efficient implementation of Edmond's
deba@327
  1961
  /// maximum weighted perfect matching algorithm. The implementation
deba@326
  1962
  /// is based on extensive use of priority queues and provides
kpeter@559
  1963
  /// \f$O(nm\log n)\f$ time complexity.
deba@326
  1964
  ///
deba@326
  1965
  /// The maximum weighted matching problem is to find undirected
deba@327
  1966
  /// edges in the graph with maximum overall weight and no two of
deba@327
  1967
  /// them shares their ends and covers all nodes. The problem can be
deba@327
  1968
  /// formulated with the following linear program.
deba@326
  1969
  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
deba@327
  1970
  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
deba@327
  1971
      \quad \forall B\in\mathcal{O}\f] */
deba@326
  1972
  /// \f[x_e \ge 0\quad \forall e\in E\f]
deba@326
  1973
  /// \f[\max \sum_{e\in E}x_ew_e\f]
deba@327
  1974
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
deba@327
  1975
  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
deba@327
  1976
  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
deba@327
  1977
  /// subsets of the nodes.
deba@326
  1978
  ///
deba@326
  1979
  /// The algorithm calculates an optimal matching and a proof of the
deba@326
  1980
  /// optimality. The solution of the dual problem can be used to check
deba@327
  1981
  /// the result of the algorithm. The dual linear problem is the
deba@327
  1982
  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
deba@327
  1983
      w_{uv} \quad \forall uv\in E\f] */
deba@326
  1984
  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
deba@327
  1985
  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
deba@327
  1986
      \frac{\vert B \vert - 1}{2}z_B\f] */
deba@326
  1987
  ///
deba@326
  1988
  /// The algorithm can be executed with \c run() or the \c init() and
deba@326
  1989
  /// then the \c start() member functions. After it the matching can
deba@326
  1990
  /// be asked with \c matching() or mate() functions. The dual
deba@326
  1991
  /// solution can be get with \c nodeValue(), \c blossomNum() and \c
deba@326
  1992
  /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
deba@326
  1993
  /// "BlossomIt" nested class which is able to iterate on the nodes
deba@326
  1994
  /// of a blossom. If the value type is integral then the dual
deba@326
  1995
  /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
kpeter@559
  1996
  template <typename GR,
kpeter@559
  1997
            typename WM = typename GR::template EdgeMap<int> >
deba@326
  1998
  class MaxWeightedPerfectMatching {
deba@326
  1999
  public:
deba@326
  2000
kpeter@559
  2001
    typedef GR Graph;
kpeter@559
  2002
    typedef WM WeightMap;
deba@326
  2003
    typedef typename WeightMap::Value Value;
deba@326
  2004
deba@326
  2005
    /// \brief Scaling factor for dual solution
deba@326
  2006
    ///
deba@326
  2007
    /// Scaling factor for dual solution, it is equal to 4 or 1
deba@326
  2008
    /// according to the value type.
deba@326
  2009
    static const int dualScale =
deba@326
  2010
      std::numeric_limits<Value>::is_integer ? 4 : 1;
deba@326
  2011
deba@326
  2012
    typedef typename Graph::template NodeMap<typename Graph::Arc>
deba@326
  2013
    MatchingMap;
deba@326
  2014
deba@326
  2015
  private:
deba@326
  2016
deba@326
  2017
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
deba@326
  2018
deba@326
  2019
    typedef typename Graph::template NodeMap<Value> NodePotential;
deba@326
  2020
    typedef std::vector<Node> BlossomNodeList;
deba@326
  2021
deba@326
  2022
    struct BlossomVariable {
deba@326
  2023
      int begin, end;
deba@326
  2024
      Value value;
deba@326
  2025
deba@326
  2026
      BlossomVariable(int _begin, int _end, Value _value)
deba@326
  2027
        : begin(_begin), end(_end), value(_value) {}
deba@326
  2028
deba@326
  2029
    };
deba@326
  2030
deba@326
  2031
    typedef std::vector<BlossomVariable> BlossomPotential;
deba@326
  2032
deba@326
  2033
    const Graph& _graph;
deba@326
  2034
    const WeightMap& _weight;
deba@326
  2035
deba@326
  2036
    MatchingMap* _matching;
deba@326
  2037
deba@326
  2038
    NodePotential* _node_potential;
deba@326
  2039
deba@326
  2040
    BlossomPotential _blossom_potential;
deba@326
  2041
    BlossomNodeList _blossom_node_list;
deba@326
  2042
deba@326
  2043
    int _node_num;
deba@326
  2044
    int _blossom_num;
deba@326
  2045
deba@326
  2046
    typedef RangeMap<int> IntIntMap;
deba@326
  2047
deba@326
  2048
    enum Status {
deba@326
  2049
      EVEN = -1, MATCHED = 0, ODD = 1
deba@326
  2050
    };
deba@326
  2051
deba@327
  2052
    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
deba@326
  2053
    struct BlossomData {
deba@326
  2054
      int tree;
deba@326
  2055
      Status status;
deba@326
  2056
      Arc pred, next;
deba@326
  2057
      Value pot, offset;
deba@326
  2058
    };
deba@326
  2059
deba@327
  2060
    IntNodeMap *_blossom_index;
deba@326
  2061
    BlossomSet *_blossom_set;
deba@326
  2062
    RangeMap<BlossomData>* _blossom_data;
deba@326
  2063
deba@327
  2064
    IntNodeMap *_node_index;
deba@327
  2065
    IntArcMap *_node_heap_index;
deba@326
  2066
deba@326
  2067
    struct NodeData {
deba@326
  2068
deba@327
  2069
      NodeData(IntArcMap& node_heap_index)
deba@326
  2070
        : heap(node_heap_index) {}
deba@326
  2071
deba@326
  2072
      int blossom;
deba@326
  2073
      Value pot;
deba@327
  2074
      BinHeap<Value, IntArcMap> heap;
deba@326
  2075
      std::map<int, Arc> heap_index;
deba@326
  2076
deba@326
  2077
      int tree;
deba@326
  2078
    };
deba@326
  2079
deba@326
  2080
    RangeMap<NodeData>* _node_data;
deba@326
  2081
deba@326
  2082
    typedef ExtendFindEnum<IntIntMap> TreeSet;
deba@326
  2083
deba@326
  2084
    IntIntMap *_tree_set_index;
deba@326
  2085
    TreeSet *_tree_set;
deba@326
  2086
deba@326
  2087
    IntIntMap *_delta2_index;
deba@326
  2088
    BinHeap<Value, IntIntMap> *_delta2;
deba@326
  2089
deba@327
  2090
    IntEdgeMap *_delta3_index;
deba@327
  2091
    BinHeap<Value, IntEdgeMap> *_delta3;
deba@326
  2092
deba@326
  2093
    IntIntMap *_delta4_index;
deba@326
  2094
    BinHeap<Value, IntIntMap> *_delta4;
deba@326
  2095
deba@326
  2096
    Value _delta_sum;
deba@326
  2097
deba@326
  2098
    void createStructures() {
deba@326
  2099
      _node_num = countNodes(_graph);
deba@326
  2100
      _blossom_num = _node_num * 3 / 2;
deba@326
  2101
deba@326
  2102
      if (!_matching) {
deba@326
  2103
        _matching = new MatchingMap(_graph);
deba@326
  2104
      }
deba@326
  2105
      if (!_node_potential) {
deba@326
  2106
        _node_potential = new NodePotential(_graph);
deba@326
  2107
      }
deba@326
  2108
      if (!_blossom_set) {
deba@327
  2109
        _blossom_index = new IntNodeMap(_graph);
deba@326
  2110
        _blossom_set = new BlossomSet(*_blossom_index);
deba@326
  2111
        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
deba@326
  2112
      }
deba@326
  2113
deba@326
  2114
      if (!_node_index) {
deba@327
  2115
        _node_index = new IntNodeMap(_graph);
deba@327
  2116
        _node_heap_index = new IntArcMap(_graph);
deba@326
  2117
        _node_data = new RangeMap<NodeData>(_node_num,
deba@327
  2118
                                            NodeData(*_node_heap_index));
deba@326
  2119
      }
deba@326
  2120
deba@326
  2121
      if (!_tree_set) {
deba@326
  2122
        _tree_set_index = new IntIntMap(_blossom_num);
deba@326
  2123
        _tree_set = new TreeSet(*_tree_set_index);
deba@326
  2124
      }
deba@326
  2125
      if (!_delta2) {
deba@326
  2126
        _delta2_index = new IntIntMap(_blossom_num);
deba@326
  2127
        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
deba@326
  2128
      }
deba@326
  2129
      if (!_delta3) {
deba@327
  2130
        _delta3_index = new IntEdgeMap(_graph);
deba@327
  2131
        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
deba@326
  2132
      }
deba@326
  2133
      if (!_delta4) {
deba@326
  2134
        _delta4_index = new IntIntMap(_blossom_num);
deba@326
  2135
        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
deba@326
  2136
      }
deba@326
  2137
    }
deba@326
  2138
deba@326
  2139
    void destroyStructures() {
deba@326
  2140
      _node_num = countNodes(_graph);
deba@326
  2141
      _blossom_num = _node_num * 3 / 2;
deba@326
  2142
deba@326
  2143
      if (_matching) {
deba@326
  2144
        delete _matching;
deba@326
  2145
      }
deba@326
  2146
      if (_node_potential) {
deba@326
  2147
        delete _node_potential;
deba@326
  2148
      }
deba@326
  2149
      if (_blossom_set) {
deba@326
  2150
        delete _blossom_index;
deba@326
  2151
        delete _blossom_set;
deba@326
  2152
        delete _blossom_data;
deba@326
  2153
      }
deba@326
  2154
deba@326
  2155
      if (_node_index) {
deba@326
  2156
        delete _node_index;
deba@326
  2157
        delete _node_heap_index;
deba@326
  2158
        delete _node_data;
deba@326
  2159
      }
deba@326
  2160
deba@326
  2161
      if (_tree_set) {
deba@326
  2162
        delete _tree_set_index;
deba@326
  2163
        delete _tree_set;
deba@326
  2164
      }
deba@326
  2165
      if (_delta2) {
deba@326
  2166
        delete _delta2_index;
deba@326
  2167
        delete _delta2;
deba@326
  2168
      }
deba@326
  2169
      if (_delta3) {
deba@326
  2170
        delete _delta3_index;
deba@326
  2171
        delete _delta3;
deba@326
  2172
      }
deba@326
  2173
      if (_delta4) {
deba@326
  2174
        delete _delta4_index;
deba@326
  2175
        delete _delta4;
deba@326
  2176
      }
deba@326
  2177
    }
deba@326
  2178
deba@326
  2179
    void matchedToEven(int blossom, int tree) {
deba@326
  2180
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
deba@326
  2181
        _delta2->erase(blossom);
deba@326
  2182
      }
deba@326
  2183
deba@326
  2184
      if (!_blossom_set->trivial(blossom)) {
deba@326
  2185
        (*_blossom_data)[blossom].pot -=
deba@326
  2186
          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
deba@326
  2187
      }
deba@326
  2188
deba@326
  2189
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
deba@326
  2190
           n != INVALID; ++n) {
deba@326
  2191
deba@326
  2192
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
deba@326
  2193
        int ni = (*_node_index)[n];
deba@326
  2194
deba@326
  2195
        (*_node_data)[ni].heap.clear();
deba@326
  2196
        (*_node_data)[ni].heap_index.clear();
deba@326
  2197
deba@326
  2198
        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
deba@326
  2199
deba@326
  2200
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
deba@326
  2201
          Node v = _graph.source(e);
deba@326
  2202
          int vb = _blossom_set->find(v);
deba@326
  2203
          int vi = (*_node_index)[v];
deba@326
  2204
deba@326
  2205
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
deba@326
  2206
            dualScale * _weight[e];
deba@326
  2207
deba@326
  2208
          if ((*_blossom_data)[vb].status == EVEN) {
deba@326
  2209
            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
deba@326
  2210
              _delta3->push(e, rw / 2);
deba@326
  2211
            }
deba@326
  2212
          } else {
deba@326
  2213
            typename std::map<int, Arc>::iterator it =
deba@326
  2214
              (*_node_data)[vi].heap_index.find(tree);
deba@326
  2215
deba@326
  2216
            if (it != (*_node_data)[vi].heap_index.end()) {
deba@326
  2217
              if ((*_node_data)[vi].heap[it->second] > rw) {
deba@326
  2218
                (*_node_data)[vi].heap.replace(it->second, e);
deba@326
  2219
                (*_node_data)[vi].heap.decrease(e, rw);
deba@326
  2220
                it->second = e;
deba@326
  2221
              }
deba@326
  2222
            } else {
deba@326
  2223
              (*_node_data)[vi].heap.push(e, rw);
deba@326
  2224
              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
deba@326
  2225
            }
deba@326
  2226
deba@326
  2227
            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
deba@326
  2228
              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
deba@326
  2229
deba@326
  2230
              if ((*_blossom_data)[vb].status == MATCHED) {
deba@326
  2231
                if (_delta2->state(vb) != _delta2->IN_HEAP) {
deba@326
  2232
                  _delta2->push(vb, _blossom_set->classPrio(vb) -
deba@326
  2233
                               (*_blossom_data)[vb].offset);
deba@326
  2234
                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
deba@326
  2235
                           (*_blossom_data)[vb].offset){
deba@326
  2236
                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
deba@326
  2237
                                   (*_blossom_data)[vb].offset);
deba@326
  2238
                }
deba@326
  2239
              }
deba@326
  2240
            }
deba@326
  2241
          }
deba@326
  2242
        }
deba@326
  2243
      }
deba@326
  2244
      (*_blossom_data)[blossom].offset = 0;
deba@326
  2245
    }
deba@326
  2246
deba@326
  2247
    void matchedToOdd(int blossom) {
deba@326
  2248
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
deba@326
  2249
        _delta2->erase(blossom);
deba@326
  2250
      }
deba@326
  2251
      (*_blossom_data)[blossom].offset += _delta_sum;
deba@326
  2252
      if (!_blossom_set->trivial(blossom)) {
deba@326
  2253
        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
deba@326
  2254
                     (*_blossom_data)[blossom].offset);
deba@326
  2255
      }
deba@326
  2256
    }
deba@326
  2257
deba@326
  2258
    void evenToMatched(int blossom, int tree) {
deba@326
  2259
      if (!_blossom_set->trivial(blossom)) {
deba@326
  2260
        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
deba@326
  2261
      }
deba@326
  2262
deba@326
  2263
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
deba@326
  2264
           n != INVALID; ++n) {
deba@326
  2265
        int ni = (*_node_index)[n];
deba@326
  2266
        (*_node_data)[ni].pot -= _delta_sum;
deba@326
  2267
deba@326
  2268
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
deba@326
  2269
          Node v = _graph.source(e);
deba@326
  2270
          int vb = _blossom_set->find(v);
deba@326
  2271
          int vi = (*_node_index)[v];
deba@326
  2272
deba@326
  2273
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
deba@326
  2274
            dualScale * _weight[e];
deba@326
  2275
deba@326
  2276
          if (vb == blossom) {
deba@326
  2277
            if (_delta3->state(e) == _delta3->IN_HEAP) {
deba@326
  2278
              _delta3->erase(e);
deba@326
  2279
            }
deba@326
  2280
          } else if ((*_blossom_data)[vb].status == EVEN) {
deba@326
  2281
deba@326
  2282
            if (_delta3->state(e) == _delta3->IN_HEAP) {
deba@326
  2283
              _delta3->erase(e);
deba@326
  2284
            }
deba@326
  2285
deba@326
  2286
            int vt = _tree_set->find(vb);
deba@326
  2287
deba@326
  2288
            if (vt != tree) {
deba@326
  2289
deba@326
  2290
              Arc r = _graph.oppositeArc(e);
deba@326
  2291
deba@326
  2292
              typename std::map<int, Arc>::iterator it =
deba@326
  2293
                (*_node_data)[ni].heap_index.find(vt);
deba@326
  2294
deba@326
  2295
              if (it != (*_node_data)[ni].heap_index.end()) {
deba@326
  2296
                if ((*_node_data)[ni].heap[it->second] > rw) {
deba@326
  2297
                  (*_node_data)[ni].heap.replace(it->second, r);
deba@326
  2298
                  (*_node_data)[ni].heap.decrease(r, rw);
deba@326
  2299
                  it->second = r;
deba@326
  2300
                }
deba@326
  2301
              } else {
deba@326
  2302
                (*_node_data)[ni].heap.push(r, rw);
deba@326
  2303
                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
deba@326
  2304
              }
deba@326
  2305
deba@326
  2306
              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
deba@326
  2307
                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
deba@326
  2308
deba@326
  2309
                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
deba@326
  2310
                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
deba@326
  2311
                               (*_blossom_data)[blossom].offset);
deba@326
  2312
                } else if ((*_delta2)[blossom] >
deba@326
  2313
                           _blossom_set->classPrio(blossom) -
deba@326
  2314
                           (*_blossom_data)[blossom].offset){
deba@326
  2315
                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
deba@326
  2316
                                   (*_blossom_data)[blossom].offset);
deba@326
  2317
                }
deba@326
  2318
              }
deba@326
  2319
            }
deba@326
  2320
          } else {
deba@326
  2321
deba@326
  2322
            typename std::map<int, Arc>::iterator it =
deba@326
  2323
              (*_node_data)[vi].heap_index.find(tree);
deba@326
  2324
deba@326
  2325
            if (it != (*_node_data)[vi].heap_index.end()) {
deba@326
  2326
              (*_node_data)[vi].heap.erase(it->second);
deba@326
  2327
              (*_node_data)[vi].heap_index.erase(it);
deba@326
  2328
              if ((*_node_data)[vi].heap.empty()) {
deba@326
  2329
                _blossom_set->increase(v, std::numeric_limits<Value>::max());
deba@326
  2330
              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
deba@326
  2331
                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
deba@326
  2332
              }
deba@326
  2333
deba@326
  2334
              if ((*_blossom_data)[vb].status == MATCHED) {
deba@326
  2335
                if (_blossom_set->classPrio(vb) ==
deba@326
  2336
                    std::numeric_limits<Value>::max()) {
deba@326
  2337
                  _delta2->erase(vb);
deba@326
  2338
                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
deba@326
  2339
                           (*_blossom_data)[vb].offset) {
deba@326
  2340
                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
deba@326
  2341
                                   (*_blossom_data)[vb].offset);
deba@326
  2342
                }
deba@326
  2343
              }
deba@326
  2344
            }
deba@326
  2345
          }
deba@326
  2346
        }
deba@326
  2347
      }
deba@326
  2348
    }
deba@326
  2349
deba@326
  2350
    void oddToMatched(int blossom) {
deba@326
  2351
      (*_blossom_data)[blossom].offset -= _delta_sum;
deba@326
  2352
deba@326
  2353
      if (_blossom_set->classPrio(blossom) !=
deba@326
  2354
          std::numeric_limits<Value>::max()) {
deba@326
  2355
        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
deba@326
  2356
                       (*_blossom_data)[blossom].offset);
deba@326
  2357
      }
deba@326
  2358
deba@326
  2359
      if (!_blossom_set->trivial(blossom)) {
deba@326
  2360
        _delta4->erase(blossom);
deba@326
  2361
      }
deba@326
  2362
    }
deba@326
  2363
deba@326
  2364
    void oddToEven(int blossom, int tree) {
deba@326
  2365
      if (!_blossom_set->trivial(blossom)) {
deba@326
  2366
        _delta4->erase(blossom);
deba@326
  2367
        (*_blossom_data)[blossom].pot -=
deba@326
  2368
          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
deba@326
  2369
      }
deba@326
  2370
deba@326
  2371
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
deba@326
  2372
           n != INVALID; ++n) {
deba@326
  2373
        int ni = (*_node_index)[n];
deba@326
  2374
deba@326
  2375
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
deba@326
  2376
deba@326
  2377
        (*_node_data)[ni].heap.clear();
deba@326
  2378
        (*_node_data)[ni].heap_index.clear();
deba@326
  2379
        (*_node_data)[ni].pot +=
deba@326
  2380
          2 * _delta_sum - (*_blossom_data)[blossom].offset;
deba@326
  2381
deba@326
  2382
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
deba@326
  2383
          Node v = _graph.source(e);
deba@326
  2384
          int vb = _blossom_set->find(v);
deba@326
  2385
          int vi = (*_node_index)[v];
deba@326
  2386
deba@326
  2387
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
deba@326
  2388
            dualScale * _weight[e];
deba@326
  2389
deba@326
  2390
          if ((*_blossom_data)[vb].status == EVEN) {
deba@326
  2391
            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
deba@326
  2392
              _delta3->push(e, rw / 2);
deba@326
  2393
            }
deba@326
  2394
          } else {
deba@326
  2395
deba@326
  2396
            typename std::map<int, Arc>::iterator it =
deba@326
  2397
              (*_node_data)[vi].heap_index.find(tree);
deba@326
  2398
deba@326
  2399
            if (it != (*_node_data)[vi].heap_index.end()) {
deba@326
  2400
              if ((*_node_data)[vi].heap[it->second] > rw) {
deba@326
  2401
                (*_node_data)[vi].heap.replace(it->second, e);
deba@326
  2402
                (*_node_data)[vi].heap.decrease(e, rw);
deba@326
  2403
                it->second = e;
deba@326
  2404
              }
deba@326
  2405
            } else {
deba@326
  2406
              (*_node_data)[vi].heap.push(e, rw);
deba@326
  2407
              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
deba@326
  2408
            }
deba@326
  2409
deba@326
  2410
            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
deba@326
  2411
              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
deba@326
  2412
deba@326
  2413
              if ((*_blossom_data)[vb].status == MATCHED) {
deba@326
  2414
                if (_delta2->state(vb) != _delta2->IN_HEAP) {
deba@326
  2415
                  _delta2->push(vb, _blossom_set->classPrio(vb) -
deba@326
  2416
                               (*_blossom_data)[vb].offset);
deba@326
  2417
                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
deba@326
  2418
                           (*_blossom_data)[vb].offset) {
deba@326
  2419
                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
deba@326
  2420
                                   (*_blossom_data)[vb].offset);
deba@326
  2421
                }
deba@326
  2422
              }
deba@326
  2423
            }
deba@326
  2424
          }
deba@326
  2425
        }
deba@326
  2426
      }
deba@326
  2427
      (*_blossom_data)[blossom].offset = 0;
deba@326
  2428
    }
deba@326
  2429
deba@326
  2430
    void alternatePath(int even, int tree) {
deba@326
  2431
      int odd;
deba@326
  2432
deba@326
  2433
      evenToMatched(even, tree);
deba@326
  2434
      (*_blossom_data)[even].status = MATCHED;
deba@326
  2435
deba@326
  2436
      while ((*_blossom_data)[even].pred != INVALID) {
deba@326
  2437
        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
deba@326
  2438
        (*_blossom_data)[odd].status = MATCHED;
deba@326
  2439
        oddToMatched(odd);
deba@326
  2440
        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
deba@326
  2441
deba@326
  2442
        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
deba@326
  2443
        (*_blossom_data)[even].status = MATCHED;
deba@326
  2444
        evenToMatched(even, tree);
deba@326
  2445
        (*_blossom_data)[even].next =
deba@326
  2446
          _graph.oppositeArc((*_blossom_data)[odd].pred);
deba@326
  2447
      }
deba@326
  2448
deba@326
  2449
    }
deba@326
  2450
deba@326
  2451
    void destroyTree(int tree) {
deba@326
  2452
      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
deba@326
  2453
        if ((*_blossom_data)[b].status == EVEN) {
deba@326
  2454
          (*_blossom_data)[b].status = MATCHED;
deba@326
  2455
          evenToMatched(b, tree);
deba@326
  2456
        } else if ((*_blossom_data)[b].status == ODD) {
deba@326
  2457
          (*_blossom_data)[b].status = MATCHED;
deba@326
  2458
          oddToMatched(b);
deba@326
  2459
        }
deba@326
  2460
      }
deba@326
  2461
      _tree_set->eraseClass(tree);
deba@326
  2462
    }
deba@326
  2463
deba@327
  2464
    void augmentOnEdge(const Edge& edge) {
deba@327
  2465
deba@327
  2466
      int left = _blossom_set->find(_graph.u(edge));
deba@327
  2467
      int right = _blossom_set->find(_graph.v(edge));
deba@326
  2468
deba@326
  2469
      int left_tree = _tree_set->find(left);
deba@326
  2470
      alternatePath(left, left_tree);
deba@326
  2471
      destroyTree(left_tree);
deba@326
  2472
deba@326
  2473
      int right_tree = _tree_set->find(right);
deba@326
  2474
      alternatePath(right, right_tree);
deba@326
  2475
      destroyTree(right_tree);
deba@326
  2476
deba@327
  2477
      (*_blossom_data)[left].next = _graph.direct(edge, true);
deba@327
  2478
      (*_blossom_data)[right].next = _graph.direct(edge, false);
deba@326
  2479
    }
deba@326
  2480
deba@326
  2481
    void extendOnArc(const Arc& arc) {
deba@326
  2482
      int base = _blossom_set->find(_graph.target(arc));
deba@326
  2483
      int tree = _tree_set->find(base);
deba@326
  2484
deba@326
  2485
      int odd = _blossom_set->find(_graph.source(arc));
deba@326
  2486
      _tree_set->insert(odd, tree);
deba@326
  2487
      (*_blossom_data)[odd].status = ODD;
deba@326
  2488
      matchedToOdd(odd);
deba@326
  2489
      (*_blossom_data)[odd].pred = arc;
deba@326
  2490
deba@326
  2491
      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
deba@326
  2492
      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
deba@326
  2493
      _tree_set->insert(even, tree);
deba@326
  2494
      (*_blossom_data)[even].status = EVEN;
deba@326
  2495
      matchedToEven(even, tree);
deba@326
  2496
    }
deba@326
  2497
deba@327
  2498
    void shrinkOnEdge(const Edge& edge, int tree) {
deba@326
  2499
      int nca = -1;
deba@326
  2500
      std::vector<int> left_path, right_path;
deba@326
  2501
deba@326
  2502
      {
deba@326
  2503
        std::set<int> left_set, right_set;
deba@326
  2504
        int left = _blossom_set->find(_graph.u(edge));
deba@326
  2505
        left_path.push_back(left);
deba@326
  2506
        left_set.insert(left);
deba@326
  2507
deba@326
  2508
        int right = _blossom_set->find(_graph.v(edge));
deba@326
  2509
        right_path.push_back(right);
deba@326
  2510
        right_set.insert(right);
deba@326
  2511
deba@326
  2512
        while (true) {
deba@326
  2513
deba@326
  2514
          if ((*_blossom_data)[left].pred == INVALID) break;
deba@326
  2515
deba@326
  2516
          left =
deba@326
  2517
            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
deba@326
  2518
          left_path.push_back(left);
deba@326
  2519
          left =
deba@326
  2520
            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
deba@326
  2521
          left_path.push_back(left);
deba@326
  2522
deba@326
  2523
          left_set.insert(left);
deba@326
  2524
deba@326
  2525
          if (right_set.find(left) != right_set.end()) {
deba@326
  2526
            nca = left;
deba@326
  2527
            break;
deba@326
  2528
          }
deba@326
  2529
deba@326
  2530
          if ((*_blossom_data)[right].pred == INVALID) break;
deba@326
  2531
deba@326
  2532
          right =
deba@326
  2533
            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
deba@326
  2534
          right_path.push_back(right);
deba@326
  2535
          right =
deba@326
  2536
            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
deba@326
  2537
          right_path.push_back(right);
deba@326
  2538
deba@326
  2539
          right_set.insert(right);
deba@326
  2540
deba@326
  2541
          if (left_set.find(right) != left_set.end()) {
deba@326
  2542
            nca = right;
deba@326
  2543
            break;
deba@326
  2544
          }
deba@326
  2545
deba@326
  2546
        }
deba@326
  2547
deba@326
  2548
        if (nca == -1) {
deba@326
  2549
          if ((*_blossom_data)[left].pred == INVALID) {
deba@326
  2550
            nca = right;
deba@326
  2551
            while (left_set.find(nca) == left_set.end()) {
deba@326
  2552
              nca =
deba@326
  2553
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
deba@326
  2554
              right_path.push_back(nca);
deba@326
  2555
              nca =
deba@326
  2556
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
deba@326
  2557
              right_path.push_back(nca);
deba@326
  2558
            }
deba@326
  2559
          } else {
deba@326
  2560
            nca = left;
deba@326
  2561
            while (right_set.find(nca) == right_set.end()) {
deba@326
  2562
              nca =
deba@326
  2563
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
deba@326
  2564
              left_path.push_back(nca);
deba@326
  2565
              nca =
deba@326
  2566
                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
deba@326
  2567
              left_path.push_back(nca);
deba@326
  2568
            }
deba@326
  2569
          }
deba@326
  2570
        }
deba@326
  2571
      }
deba@326
  2572
deba@326
  2573
      std::vector<int> subblossoms;
deba@326
  2574
      Arc prev;
deba@326
  2575
deba@326
  2576
      prev = _graph.direct(edge, true);
deba@326
  2577
      for (int i = 0; left_path[i] != nca; i += 2) {
deba@326
  2578
        subblossoms.push_back(left_path[i]);
deba@326
  2579
        (*_blossom_data)[left_path[i]].next = prev;
deba@326
  2580
        _tree_set->erase(left_path[i]);
deba@326
  2581
deba@326
  2582
        subblossoms.push_back(left_path[i + 1]);
deba@326
  2583
        (*_blossom_data)[left_path[i + 1]].status = EVEN;
deba@326
  2584
        oddToEven(left_path[i + 1], tree);
deba@326
  2585
        _tree_set->erase(left_path[i + 1]);
deba@326
  2586
        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
deba@326
  2587
      }
deba@326
  2588
deba@326
  2589
      int k = 0;
deba@326
  2590
      while (right_path[k] != nca) ++k;
deba@326
  2591
deba@326
  2592
      subblossoms.push_back(nca);
deba@326
  2593
      (*_blossom_data)[nca].next = prev;
deba@326
  2594
deba@326
  2595
      for (int i = k - 2; i >= 0; i -= 2) {
deba@326
  2596
        subblossoms.push_back(right_path[i + 1]);
deba@326
  2597
        (*_blossom_data)[right_path[i + 1]].status = EVEN;
deba@326
  2598
        oddToEven(right_path[i + 1], tree);
deba@326
  2599
        _tree_set->erase(right_path[i + 1]);
deba@326
  2600
deba@326
  2601
        (*_blossom_data)[right_path[i + 1]].next =
deba@326
  2602
          (*_blossom_data)[right_path[i + 1]].pred;
deba@326
  2603
deba@326
  2604
        subblossoms.push_back(right_path[i]);
deba@326
  2605
        _tree_set->erase(right_path[i]);
deba@326
  2606
      }
deba@326
  2607
deba@326
  2608
      int surface =
deba@326
  2609
        _blossom_set->join(subblossoms.begin(), subblossoms.end());
deba@326
  2610
deba@326
  2611
      for (int i = 0; i < int(subblossoms.size()); ++i) {
deba@326
  2612
        if (!_blossom_set->trivial(subblossoms[i])) {
deba@326
  2613
          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
deba@326
  2614
        }
deba@326
  2615
        (*_blossom_data)[subblossoms[i]].status = MATCHED;
deba@326
  2616
      }
deba@326
  2617
deba@326
  2618
      (*_blossom_data)[surface].pot = -2 * _delta_sum;
deba@326
  2619
      (*_blossom_data)[surface].offset = 0;
deba@326
  2620
      (*_blossom_data)[surface].status = EVEN;
deba@326
  2621
      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
deba@326
  2622
      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
deba@326
  2623
deba@326
  2624
      _tree_set->insert(surface, tree);
deba@326
  2625
      _tree_set->erase(nca);
deba@326
  2626
    }
deba@326
  2627
deba@326
  2628
    void splitBlossom(int blossom) {
deba@326
  2629
      Arc next = (*_blossom_data)[blossom].next;
deba@326
  2630
      Arc pred = (*_blossom_data)[blossom].pred;
deba@326
  2631
deba@326
  2632
      int tree = _tree_set->find(blossom);
deba@326
  2633
deba@326
  2634
      (*_blossom_data)[blossom].status = MATCHED;
deba@326
  2635
      oddToMatched(blossom);
deba@326
  2636
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
deba@326
  2637
        _delta2->erase(blossom);
deba@326
  2638
      }
deba@326
  2639
deba@326
  2640
      std::vector<int> subblossoms;
deba@326
  2641
      _blossom_set->split(blossom, std::back_inserter(subblossoms));
deba@326
  2642
deba@326
  2643
      Value offset = (*_blossom_data)[blossom].offset;
deba@326
  2644
      int b = _blossom_set->find(_graph.source(pred));
deba@326
  2645
      int d = _blossom_set->find(_graph.source(next));
deba@326
  2646
deba@326
  2647
      int ib = -1, id = -1;
deba@326
  2648
      for (int i = 0; i < int(subblossoms.size()); ++i) {
deba@326
  2649
        if (subblossoms[i] == b) ib = i;
deba@326
  2650
        if (subblossoms[i] == d) id = i;
deba@326
  2651
deba@326
  2652
        (*_blossom_data)[subblossoms[i]].offset = offset;
deba@326
  2653
        if (!_blossom_set->trivial(subblossoms[i])) {
deba@326
  2654
          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
deba@326
  2655
        }
deba@326
  2656
        if (_blossom_set->classPrio(subblossoms[i]) !=
deba@326
  2657
            std::numeric_limits<Value>::max()) {
deba@326
  2658
          _delta2->push(subblossoms[i],
deba@326
  2659
                        _blossom_set->classPrio(subblossoms[i]) -
deba@326
  2660
                        (*_blossom_data)[subblossoms[i]].offset);
deba@326
  2661
        }
deba@326
  2662
      }
deba@326
  2663
deba@326
  2664
      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
deba@326
  2665
        for (int i = (id + 1) % subblossoms.size();
deba@326
  2666
             i != ib; i = (i + 2) % subblossoms.size()) {
deba@326
  2667
          int sb = subblossoms[i];
deba@326
  2668
          int tb = subblossoms[(i + 1) % subblossoms.size()];
deba@326
  2669
          (*_blossom_data)[sb].next =
deba@326
  2670
            _graph.oppositeArc((*_blossom_data)[tb].next);
deba@326
  2671
        }
deba@326
  2672
deba@326
  2673
        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
deba@326
  2674
          int sb = subblossoms[i];
deba@326
  2675
          int tb = subblossoms[(i + 1) % subblossoms.size()];
deba@326
  2676
          int ub = subblossoms[(i + 2) % subblossoms.size()];
deba@326
  2677
deba@326
  2678
          (*_blossom_data)[sb].status = ODD;
deba@326
  2679
          matchedToOdd(sb);
deba@326
  2680
          _tree_set->insert(sb, tree);
deba@326
  2681
          (*_blossom_data)[sb].pred = pred;
deba@326
  2682
          (*_blossom_data)[sb].next =
deba@326
  2683
                           _graph.oppositeArc((*_blossom_data)[tb].next);
deba@326
  2684
deba@326
  2685
          pred = (*_blossom_data)[ub].next;
deba@326
  2686
deba@326
  2687
          (*_blossom_data)[tb].status = EVEN;
deba@326
  2688
          matchedToEven(tb, tree);
deba@326
  2689
          _tree_set->insert(tb, tree);
deba@326
  2690
          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
deba@326
  2691
        }
deba@326
  2692
deba@326
  2693
        (*_blossom_data)[subblossoms[id]].status = ODD;
deba@326
  2694
        matchedToOdd(subblossoms[id]);
deba@326
  2695
        _tree_set->insert(subblossoms[id], tree);
deba@326
  2696
        (*_blossom_data)[subblossoms[id]].next = next;
deba@326
  2697
        (*_blossom_data)[subblossoms[id]].pred = pred;
deba@326
  2698
deba@326
  2699
      } else {
deba@326
  2700
deba@326
  2701
        for (int i = (ib + 1) % subblossoms.size();
deba@326
  2702
             i != id; i = (i + 2) % subblossoms.size()) {
deba@326
  2703
          int sb = subblossoms[i];
deba@326
  2704
          int tb = subblossoms[(i + 1) % subblossoms.size()];
deba@326
  2705
          (*_blossom_data)[sb].next =
deba@326
  2706
            _graph.oppositeArc((*_blossom_data)[tb].next);
deba@326
  2707
        }
deba@326
  2708
deba@326
  2709
        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
deba@326
  2710
          int sb = subblossoms[i];
deba@326
  2711
          int tb = subblossoms[(i + 1) % subblossoms.size()];
deba@326
  2712
          int ub = subblossoms[(i + 2) % subblossoms.size()];
deba@326
  2713
deba@326
  2714
          (*_blossom_data)[sb].status = ODD;
deba@326
  2715
          matchedToOdd(sb);
deba@326
  2716
          _tree_set->insert(sb, tree);
deba@326
  2717
          (*_blossom_data)[sb].next = next;
deba@326
  2718
          (*_blossom_data)[sb].pred =
deba@326
  2719
            _graph.oppositeArc((*_blossom_data)[tb].next);
deba@326
  2720
deba@326
  2721
          (*_blossom_data)[tb].status = EVEN;
deba@326
  2722
          matchedToEven(tb, tree);
deba@326
  2723
          _tree_set->insert(tb, tree);
deba@326
  2724
          (*_blossom_data)[tb].pred =
deba@326
  2725
            (*_blossom_data)[tb].next =
deba@326
  2726
            _graph.oppositeArc((*_blossom_data)[ub].next);
deba@326
  2727
          next = (*_blossom_data)[ub].next;
deba@326
  2728
        }
deba@326
  2729
deba@326
  2730
        (*_blossom_data)[subblossoms[ib]].status = ODD;
deba@326
  2731
        matchedToOdd(subblossoms[ib]);
deba@326
  2732
        _tree_set->insert(subblossoms[ib], tree);
deba@326
  2733
        (*_blossom_data)[subblossoms[ib]].next = next;
deba@326
  2734
        (*_blossom_data)[subblossoms[ib]].pred = pred;
deba@326
  2735
      }
deba@326
  2736
      _tree_set->erase(blossom);
deba@326
  2737
    }
deba@326
  2738
deba@326
  2739
    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
deba@326
  2740
      if (_blossom_set->trivial(blossom)) {
deba@326
  2741
        int bi = (*_node_index)[base];
deba@326
  2742
        Value pot = (*_node_data)[bi].pot;
deba@326
  2743
kpeter@581
  2744
        (*_matching)[base] = matching;
deba@326
  2745
        _blossom_node_list.push_back(base);
kpeter@581
  2746
        (*_node_potential)[base] = pot;
deba@326
  2747
      } else {
deba@326
  2748
deba@326
  2749
        Value pot = (*_blossom_data)[blossom].pot;
deba@326
  2750
        int bn = _blossom_node_list.size();
deba@326
  2751
deba@326
  2752
        std::vector<int> subblossoms;
deba@326
  2753
        _blossom_set->split(blossom, std::back_inserter(subblossoms));
deba@326
  2754
        int b = _blossom_set->find(base);
deba@326
  2755
        int ib = -1;
deba@326
  2756
        for (int i = 0; i < int(subblossoms.size()); ++i) {
deba@326
  2757
          if (subblossoms[i] == b) { ib = i; break; }
deba@326
  2758
        }
deba@326
  2759
deba@326
  2760
        for (int i = 1; i < int(subblossoms.size()); i += 2) {
deba@326
  2761
          int sb = subblossoms[(ib + i) % subblossoms.size()];
deba@326
  2762
          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
deba@326
  2763
deba@326
  2764
          Arc m = (*_blossom_data)[tb].next;
deba@326
  2765
          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
deba@326
  2766
          extractBlossom(tb, _graph.source(m), m);
deba@326
  2767
        }
deba@326
  2768
        extractBlossom(subblossoms[ib], base, matching);
deba@326
  2769
deba@326
  2770
        int en = _blossom_node_list.size();
deba@326
  2771
deba@326
  2772
        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
deba@326
  2773
      }
deba@326
  2774
    }
deba@326
  2775
deba@326
  2776
    void extractMatching() {
deba@326
  2777
      std::vector<int> blossoms;
deba@326
  2778
      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
deba@326
  2779
        blossoms.push_back(c);
deba@326
  2780
      }
deba@326
  2781
deba@326
  2782
      for (int i = 0; i < int(blossoms.size()); ++i) {
deba@326
  2783
deba@326
  2784
        Value offset = (*_blossom_data)[blossoms[i]].offset;
deba@326
  2785
        (*_blossom_data)[blossoms[i]].pot += 2 * offset;
deba@326
  2786
        for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
deba@326
  2787
             n != INVALID; ++n) {
deba@326
  2788
          (*_node_data)[(*_node_index)[n]].pot -= offset;
deba@326
  2789
        }
deba@326
  2790
deba@326
  2791
        Arc matching = (*_blossom_data)[blossoms[i]].next;
deba@326
  2792
        Node base = _graph.source(matching);
deba@326
  2793
        extractBlossom(blossoms[i], base, matching);
deba@326
  2794
      }
deba@326
  2795
    }
deba@326
  2796
deba@326
  2797
  public:
deba@326
  2798
deba@326
  2799
    /// \brief Constructor
deba@326
  2800
    ///
deba@326
  2801
    /// Constructor.
deba@326
  2802
    MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
deba@326
  2803
      : _graph(graph), _weight(weight), _matching(0),
deba@326
  2804
        _node_potential(0), _blossom_potential(), _blossom_node_list(),
deba@326
  2805
        _node_num(0), _blossom_num(0),
deba@326
  2806
deba@326
  2807
        _blossom_index(0), _blossom_set(0), _blossom_data(0),
deba@326
  2808
        _node_index(0), _node_heap_index(0), _node_data(0),
deba@326
  2809
        _tree_set_index(0), _tree_set(0),
deba@326
  2810
deba@326
  2811
        _delta2_index(0), _delta2(0),
deba@326
  2812
        _delta3_index(0), _delta3(0),
deba@326
  2813
        _delta4_index(0), _delta4(0),
deba@326
  2814
deba@326
  2815
        _delta_sum() {}
deba@326
  2816
deba@326
  2817
    ~MaxWeightedPerfectMatching() {
deba@326
  2818
      destroyStructures();
deba@326
  2819
    }
deba@326
  2820
deba@326
  2821
    /// \name Execution control
alpar@330
  2822
    /// The simplest way to execute the algorithm is to use the
deba@326
  2823
    /// \c run() member function.
deba@326
  2824
deba@326
  2825
    ///@{
deba@326
  2826
deba@326
  2827
    /// \brief Initialize the algorithm
deba@326
  2828
    ///
deba@326
  2829
    /// Initialize the algorithm
deba@326
  2830
    void init() {
deba@326
  2831
      createStructures();
deba@326
  2832
deba@326
  2833
      for (ArcIt e(_graph); e != INVALID; ++e) {
kpeter@581
  2834
        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
deba@326
  2835
      }
deba@326
  2836
      for (EdgeIt e(_graph); e != INVALID; ++e) {
kpeter@581
  2837
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
deba@326
  2838
      }
deba@326
  2839
      for (int i = 0; i < _blossom_num; ++i) {
kpeter@581
  2840
        (*_delta2_index)[i] = _delta2->PRE_HEAP;
kpeter@581
  2841
        (*_delta4_index)[i] = _delta4->PRE_HEAP;
deba@326
  2842
      }
deba@326
  2843
deba@326
  2844
      int index = 0;
deba@326
  2845
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@326
  2846
        Value max = - std::numeric_limits<Value>::max();
deba@326
  2847
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
deba@326
  2848
          if (_graph.target(e) == n) continue;
deba@326
  2849
          if ((dualScale * _weight[e]) / 2 > max) {
deba@326
  2850
            max = (dualScale * _weight[e]) / 2;
deba@326
  2851
          }
deba@326
  2852
        }
kpeter@581
  2853
        (*_node_index)[n] = index;
deba@326
  2854
        (*_node_data)[index].pot = max;
deba@326
  2855
        int blossom =
deba@326
  2856
          _blossom_set->insert(n, std::numeric_limits<Value>::max());
deba@326
  2857
deba@326
  2858
        _tree_set->insert(blossom);
deba@326
  2859
deba@326
  2860
        (*_blossom_data)[blossom].status = EVEN;
deba@326
  2861
        (*_blossom_data)[blossom].pred = INVALID;
deba@326
  2862
        (*_blossom_data)[blossom].next = INVALID;
deba@326
  2863
        (*_blossom_data)[blossom].pot = 0;
deba@326
  2864
        (*_blossom_data)[blossom].offset = 0;
deba@326
  2865
        ++index;
deba@326
  2866
      }
deba@326
  2867
      for (EdgeIt e(_graph); e != INVALID; ++e) {
deba@326
  2868
        int si = (*_node_index)[_graph.u(e)];
deba@326
  2869
        int ti = (*_node_index)[_graph.v(e)];
deba@326
  2870
        if (_graph.u(e) != _graph.v(e)) {
deba@326
  2871
          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
deba@326
  2872
                            dualScale * _weight[e]) / 2);
deba@326
  2873
        }
deba@326
  2874
      }
deba@326
  2875
    }
deba@326
  2876
deba@326
  2877
    /// \brief Starts the algorithm
deba@326
  2878
    ///
deba@326
  2879
    /// Starts the algorithm
deba@326
  2880
    bool start() {
deba@326
  2881
      enum OpType {
deba@326
  2882
        D2, D3, D4
deba@326
  2883
      };
deba@326
  2884
deba@326
  2885
      int unmatched = _node_num;
deba@326
  2886
      while (unmatched > 0) {
deba@326
  2887
        Value d2 = !_delta2->empty() ?
deba@326
  2888
          _delta2->prio() : std::numeric_limits<Value>::max();
deba@326
  2889
deba@326
  2890
        Value d3 = !_delta3->empty() ?
deba@326
  2891
          _delta3->prio() : std::numeric_limits<Value>::max();
deba@326
  2892
deba@326
  2893
        Value d4 = !_delta4->empty() ?
deba@326
  2894
          _delta4->prio() : std::numeric_limits<Value>::max();
deba@326
  2895
deba@326
  2896
        _delta_sum = d2; OpType ot = D2;
deba@326
  2897
        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
deba@326
  2898
        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
deba@326
  2899
deba@326
  2900
        if (_delta_sum == std::numeric_limits<Value>::max()) {
deba@326
  2901
          return false;
deba@326
  2902
        }
deba@326
  2903
deba@326
  2904
        switch (ot) {
deba@326
  2905
        case D2:
deba@326
  2906
          {
deba@326
  2907
            int blossom = _delta2->top();
deba@326
  2908
            Node n = _blossom_set->classTop(blossom);
deba@326
  2909
            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
deba@326
  2910
            extendOnArc(e);
deba@326
  2911
          }
deba@326
  2912
          break;
deba@326
  2913
        case D3:
deba@326
  2914
          {
deba@326
  2915
            Edge e = _delta3->top();
deba@326
  2916
deba@326
  2917
            int left_blossom = _blossom_set->find(_graph.u(e));
deba@326
  2918
            int right_blossom = _blossom_set->find(_graph.v(e));
deba@326
  2919
deba@326
  2920
            if (left_blossom == right_blossom) {
deba@326
  2921
              _delta3->pop();
deba@326
  2922
            } else {
deba@326
  2923
              int left_tree = _tree_set->find(left_blossom);
deba@326
  2924
              int right_tree = _tree_set->find(right_blossom);
deba@326
  2925
deba@326
  2926
              if (left_tree == right_tree) {
deba@327
  2927
                shrinkOnEdge(e, left_tree);
deba@326
  2928
              } else {
deba@327
  2929
                augmentOnEdge(e);
deba@326
  2930
                unmatched -= 2;
deba@326
  2931
              }
deba@326
  2932
            }
deba@326
  2933
          } break;
deba@326
  2934
        case D4:
deba@326
  2935
          splitBlossom(_delta4->top());
deba@326
  2936
          break;
deba@326
  2937
        }
deba@326
  2938
      }
deba@326
  2939
      extractMatching();
deba@326
  2940
      return true;
deba@326
  2941
    }
deba@326
  2942
deba@326
  2943
    /// \brief Runs %MaxWeightedPerfectMatching algorithm.
deba@326
  2944
    ///
deba@326
  2945
    /// This method runs the %MaxWeightedPerfectMatching algorithm.
deba@326
  2946
    ///
deba@326
  2947
    /// \note mwm.run() is just a shortcut of the following code.
deba@326
  2948
    /// \code
deba@326
  2949
    ///   mwm.init();
deba@326
  2950
    ///   mwm.start();
deba@326
  2951
    /// \endcode
deba@326
  2952
    bool run() {
deba@326
  2953
      init();
deba@326
  2954
      return start();
deba@326
  2955
    }
deba@326
  2956
deba@326
  2957
    /// @}
deba@326
  2958
deba@326
  2959
    /// \name Primal solution
alpar@330
  2960
    /// Functions to get the primal solution, ie. the matching.
deba@326
  2961
deba@326
  2962
    /// @{
deba@326
  2963
deba@326
  2964
    /// \brief Returns the matching value.
deba@326
  2965
    ///
deba@326
  2966
    /// Returns the matching value.
deba@326
  2967
    Value matchingValue() const {
deba@326
  2968
      Value sum = 0;
deba@326
  2969
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@326
  2970
        if ((*_matching)[n] != INVALID) {
deba@326
  2971
          sum += _weight[(*_matching)[n]];
deba@326
  2972
        }
deba@326
  2973
      }
deba@326
  2974
      return sum /= 2;
deba@326
  2975
    }
deba@326
  2976
deba@327
  2977
    /// \brief Returns true when the edge is in the matching.
deba@326
  2978
    ///
deba@327
  2979
    /// Returns true when the edge is in the matching.
deba@327
  2980
    bool matching(const Edge& edge) const {
deba@327
  2981
      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
deba@326
  2982
    }
deba@326
  2983
deba@327
  2984
    /// \brief Returns the incident matching edge.
deba@326
  2985
    ///
deba@327
  2986
    /// Returns the incident matching arc from given edge.
deba@326
  2987
    Arc matching(const Node& node) const {
deba@326
  2988
      return (*_matching)[node];
deba@326
  2989
    }
deba@326
  2990
deba@326
  2991
    /// \brief Returns the mate of the node.
deba@326
  2992
    ///
deba@326
  2993
    /// Returns the adjancent node in a mathcing arc.
deba@326
  2994
    Node mate(const Node& node) const {
deba@326
  2995
      return _graph.target((*_matching)[node]);
deba@326
  2996
    }
deba@326
  2997
deba@326
  2998
    /// @}
deba@326
  2999
deba@326
  3000
    /// \name Dual solution
alpar@330
  3001
    /// Functions to get the dual solution.
deba@326
  3002
deba@326
  3003
    /// @{
deba@326
  3004
deba@326
  3005
    /// \brief Returns the value of the dual solution.
deba@326
  3006
    ///
deba@326
  3007
    /// Returns the value of the dual solution. It should be equal to
deba@326
  3008
    /// the primal value scaled by \ref dualScale "dual scale".
deba@326
  3009
    Value dualValue() const {
deba@326
  3010
      Value sum = 0;
deba@326
  3011
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@326
  3012
        sum += nodeValue(n);
deba@326
  3013
      }
deba@326
  3014
      for (int i = 0; i < blossomNum(); ++i) {
deba@326
  3015
        sum += blossomValue(i) * (blossomSize(i) / 2);
deba@326
  3016
      }
deba@326
  3017
      return sum;
deba@326
  3018
    }
deba@326
  3019
deba@326
  3020
    /// \brief Returns the value of the node.
deba@326
  3021
    ///
deba@326
  3022
    /// Returns the the value of the node.
deba@326
  3023
    Value nodeValue(const Node& n) const {
deba@326
  3024
      return (*_node_potential)[n];
deba@326
  3025
    }
deba@326
  3026
deba@326
  3027
    /// \brief Returns the number of the blossoms in the basis.
deba@326
  3028
    ///
deba@326
  3029
    /// Returns the number of the blossoms in the basis.
deba@326
  3030
    /// \see BlossomIt
deba@326
  3031
    int blossomNum() const {
deba@326
  3032
      return _blossom_potential.size();
deba@326
  3033
    }
deba@326
  3034
deba@326
  3035
deba@326
  3036
    /// \brief Returns the number of the nodes in the blossom.
deba@326
  3037
    ///
deba@326
  3038
    /// Returns the number of the nodes in the blossom.
deba@326
  3039
    int blossomSize(int k) const {
deba@326
  3040
      return _blossom_potential[k].end - _blossom_potential[k].begin;
deba@326
  3041
    }
deba@326
  3042
deba@326
  3043
    /// \brief Returns the value of the blossom.
deba@326
  3044
    ///
deba@326
  3045
    /// Returns the the value of the blossom.
deba@326
  3046
    /// \see BlossomIt
deba@326
  3047
    Value blossomValue(int k) const {
deba@326
  3048
      return _blossom_potential[k].value;
deba@326
  3049
    }
deba@326
  3050
alpar@330
  3051
    /// \brief Iterator for obtaining the nodes of the blossom.
deba@326
  3052
    ///
alpar@330
  3053
    /// Iterator for obtaining the nodes of the blossom. This class
alpar@330
  3054
    /// provides a common lemon style iterator for listing a
deba@326
  3055
    /// subset of the nodes.
deba@326
  3056
    class BlossomIt {
deba@326
  3057
    public:
deba@326
  3058
deba@326
  3059
      /// \brief Constructor.
deba@326
  3060
      ///
alpar@330
  3061
      /// Constructor to get the nodes of the variable.
deba@326
  3062
      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
deba@326
  3063
        : _algorithm(&algorithm)
deba@326
  3064
      {
deba@326
  3065
        _index = _algorithm->_blossom_potential[variable].begin;
deba@326
  3066
        _last = _algorithm->_blossom_potential[variable].end;
deba@326
  3067
      }
deba@326
  3068
deba@326
  3069
      /// \brief Conversion to node.
deba@326
  3070
      ///
deba@326
  3071
      /// Conversion to node.
deba@326
  3072
      operator Node() const {
deba@327
  3073
        return _algorithm->_blossom_node_list[_index];
deba@326
  3074
      }
deba@326
  3075
deba@326
  3076
      /// \brief Increment operator.
deba@326
  3077
      ///
deba@326
  3078
      /// Increment operator.
deba@326
  3079
      BlossomIt& operator++() {
deba@326
  3080
        ++_index;
deba@326
  3081
        return *this;
deba@326
  3082
      }
deba@326
  3083
deba@327
  3084
      /// \brief Validity checking
deba@327
  3085
      ///
deba@327
  3086
      /// Checks whether the iterator is invalid.
deba@327
  3087
      bool operator==(Invalid) const { return _index == _last; }
deba@327
  3088
deba@327
  3089
      /// \brief Validity checking
deba@327
  3090
      ///
deba@327
  3091
      /// Checks whether the iterator is valid.
deba@327
  3092
      bool operator!=(Invalid) const { return _index != _last; }
deba@326
  3093
deba@326
  3094
    private:
deba@326
  3095
      const MaxWeightedPerfectMatching* _algorithm;
deba@326
  3096
      int _last;
deba@326
  3097
      int _index;
deba@326
  3098
    };
deba@326
  3099
deba@326
  3100
    /// @}
deba@326
  3101
deba@326
  3102
  };
deba@326
  3103
deba@326
  3104
deba@326
  3105
} //END OF NAMESPACE LEMON
deba@326
  3106
deba@326
  3107
#endif //LEMON_MAX_MATCHING_H