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