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