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