lemon/fractional_matching.h
author Balazs Dezso <deba@inf.elte.hu>
Thu, 04 Mar 2010 10:17:02 +0100
changeset 950 86613aa28a0c
parent 948 636dadefe1e6
child 951 41d7ac528c3a
permissions -rw-r--r--
Fix documentation issues (#314)
deba@948
     1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
deba@948
     2
 *
deba@948
     3
 * This file is a part of LEMON, a generic C++ optimization library.
deba@948
     4
 *
deba@948
     5
 * Copyright (C) 2003-2009
deba@948
     6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
deba@948
     7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
deba@948
     8
 *
deba@948
     9
 * Permission to use, modify and distribute this software is granted
deba@948
    10
 * provided that this copyright notice appears in all copies. For
deba@948
    11
 * precise terms see the accompanying LICENSE file.
deba@948
    12
 *
deba@948
    13
 * This software is provided "AS IS" with no warranty of any kind,
deba@948
    14
 * express or implied, and with no claim as to its suitability for any
deba@948
    15
 * purpose.
deba@948
    16
 *
deba@948
    17
 */
deba@948
    18
deba@948
    19
#ifndef LEMON_FRACTIONAL_MATCHING_H
deba@948
    20
#define LEMON_FRACTIONAL_MATCHING_H
deba@948
    21
deba@948
    22
#include <vector>
deba@948
    23
#include <queue>
deba@948
    24
#include <set>
deba@948
    25
#include <limits>
deba@948
    26
deba@948
    27
#include <lemon/core.h>
deba@948
    28
#include <lemon/unionfind.h>
deba@948
    29
#include <lemon/bin_heap.h>
deba@948
    30
#include <lemon/maps.h>
deba@948
    31
#include <lemon/assert.h>
deba@948
    32
#include <lemon/elevator.h>
deba@948
    33
deba@948
    34
///\ingroup matching
deba@948
    35
///\file
deba@948
    36
///\brief Fractional matching algorithms in general graphs.
deba@948
    37
deba@948
    38
namespace lemon {
deba@948
    39
deba@948
    40
  /// \brief Default traits class of MaxFractionalMatching class.
deba@948
    41
  ///
deba@948
    42
  /// Default traits class of MaxFractionalMatching class.
deba@948
    43
  /// \tparam GR Graph type.
deba@948
    44
  template <typename GR>
deba@948
    45
  struct MaxFractionalMatchingDefaultTraits {
deba@948
    46
deba@948
    47
    /// \brief The type of the graph the algorithm runs on.
deba@948
    48
    typedef GR Graph;
deba@948
    49
deba@948
    50
    /// \brief The type of the map that stores the matching.
deba@948
    51
    ///
deba@948
    52
    /// The type of the map that stores the matching arcs.
deba@948
    53
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
deba@948
    54
    typedef typename Graph::template NodeMap<typename GR::Arc> MatchingMap;
deba@948
    55
deba@948
    56
    /// \brief Instantiates a MatchingMap.
deba@948
    57
    ///
deba@948
    58
    /// This function instantiates a \ref MatchingMap.
deba@948
    59
    /// \param graph The graph for which we would like to define
deba@948
    60
    /// the matching map.
deba@948
    61
    static MatchingMap* createMatchingMap(const Graph& graph) {
deba@948
    62
      return new MatchingMap(graph);
deba@948
    63
    }
deba@948
    64
deba@948
    65
    /// \brief The elevator type used by MaxFractionalMatching algorithm.
deba@948
    66
    ///
deba@948
    67
    /// The elevator type used by MaxFractionalMatching algorithm.
deba@948
    68
    ///
deba@948
    69
    /// \sa Elevator
deba@948
    70
    /// \sa LinkedElevator
deba@948
    71
    typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
deba@948
    72
deba@948
    73
    /// \brief Instantiates an Elevator.
deba@948
    74
    ///
deba@948
    75
    /// This function instantiates an \ref Elevator.
deba@948
    76
    /// \param graph The graph for which we would like to define
deba@948
    77
    /// the elevator.
deba@948
    78
    /// \param max_level The maximum level of the elevator.
deba@948
    79
    static Elevator* createElevator(const Graph& graph, int max_level) {
deba@948
    80
      return new Elevator(graph, max_level);
deba@948
    81
    }
deba@948
    82
  };
deba@948
    83
deba@948
    84
  /// \ingroup matching
deba@948
    85
  ///
deba@948
    86
  /// \brief Max cardinality fractional matching
deba@948
    87
  ///
deba@948
    88
  /// This class provides an implementation of fractional matching
deba@948
    89
  /// algorithm based on push-relabel principle.
deba@948
    90
  ///
deba@948
    91
  /// The maximum cardinality fractional matching is a relaxation of the
deba@948
    92
  /// maximum cardinality matching problem where the odd set constraints
deba@948
    93
  /// are omitted.
deba@948
    94
  /// It can be formulated with the following linear program.
deba@948
    95
  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
deba@948
    96
  /// \f[x_e \ge 0\quad \forall e\in E\f]
deba@948
    97
  /// \f[\max \sum_{e\in E}x_e\f]
deba@948
    98
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
deba@948
    99
  /// \f$X\f$. The result can be represented as the union of a
deba@948
   100
  /// matching with one value edges and a set of odd length cycles
deba@948
   101
  /// with half value edges.
deba@948
   102
  ///
deba@948
   103
  /// The algorithm calculates an optimal fractional matching and a
deba@948
   104
  /// barrier. The number of adjacents of any node set minus the size
deba@948
   105
  /// of node set is a lower bound on the uncovered nodes in the
deba@948
   106
  /// graph. For maximum matching a barrier is computed which
deba@948
   107
  /// maximizes this difference.
deba@948
   108
  ///
deba@948
   109
  /// The algorithm can be executed with the run() function.  After it
deba@948
   110
  /// the matching (the primal solution) and the barrier (the dual
deba@948
   111
  /// solution) can be obtained using the query functions.
deba@948
   112
  ///
deba@948
   113
  /// The primal solution is multiplied by
deba@950
   114
  /// \ref MaxFractionalMatching::primalScale "2".
deba@948
   115
  ///
deba@948
   116
  /// \tparam GR The undirected graph type the algorithm runs on.
deba@948
   117
#ifdef DOXYGEN
deba@948
   118
  template <typename GR, typename TR>
deba@948
   119
#else
deba@948
   120
  template <typename GR,
deba@948
   121
            typename TR = MaxFractionalMatchingDefaultTraits<GR> >
deba@948
   122
#endif
deba@948
   123
  class MaxFractionalMatching {
deba@948
   124
  public:
deba@948
   125
deba@948
   126
    /// \brief The \ref MaxFractionalMatchingDefaultTraits "traits
deba@948
   127
    /// class" of the algorithm.
deba@948
   128
    typedef TR Traits;
deba@948
   129
    /// The type of the graph the algorithm runs on.
deba@948
   130
    typedef typename TR::Graph Graph;
deba@948
   131
    /// The type of the matching map.
deba@948
   132
    typedef typename TR::MatchingMap MatchingMap;
deba@948
   133
    /// The type of the elevator.
deba@948
   134
    typedef typename TR::Elevator Elevator;
deba@948
   135
deba@948
   136
    /// \brief Scaling factor for primal solution
deba@948
   137
    ///
deba@948
   138
    /// Scaling factor for primal solution.
deba@948
   139
    static const int primalScale = 2;
deba@948
   140
deba@948
   141
  private:
deba@948
   142
deba@948
   143
    const Graph &_graph;
deba@948
   144
    int _node_num;
deba@948
   145
    bool _allow_loops;
deba@948
   146
    int _empty_level;
deba@948
   147
deba@948
   148
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
deba@948
   149
deba@948
   150
    bool _local_matching;
deba@948
   151
    MatchingMap *_matching;
deba@948
   152
deba@948
   153
    bool _local_level;
deba@948
   154
    Elevator *_level;
deba@948
   155
deba@948
   156
    typedef typename Graph::template NodeMap<int> InDegMap;
deba@948
   157
    InDegMap *_indeg;
deba@948
   158
deba@948
   159
    void createStructures() {
deba@948
   160
      _node_num = countNodes(_graph);
deba@948
   161
deba@948
   162
      if (!_matching) {
deba@948
   163
        _local_matching = true;
deba@948
   164
        _matching = Traits::createMatchingMap(_graph);
deba@948
   165
      }
deba@948
   166
      if (!_level) {
deba@948
   167
        _local_level = true;
deba@948
   168
        _level = Traits::createElevator(_graph, _node_num);
deba@948
   169
      }
deba@948
   170
      if (!_indeg) {
deba@948
   171
        _indeg = new InDegMap(_graph);
deba@948
   172
      }
deba@948
   173
    }
deba@948
   174
deba@948
   175
    void destroyStructures() {
deba@948
   176
      if (_local_matching) {
deba@948
   177
        delete _matching;
deba@948
   178
      }
deba@948
   179
      if (_local_level) {
deba@948
   180
        delete _level;
deba@948
   181
      }
deba@948
   182
      if (_indeg) {
deba@948
   183
        delete _indeg;
deba@948
   184
      }
deba@948
   185
    }
deba@948
   186
deba@948
   187
    void postprocessing() {
deba@948
   188
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
   189
        if ((*_indeg)[n] != 0) continue;
deba@948
   190
        _indeg->set(n, -1);
deba@948
   191
        Node u = n;
deba@948
   192
        while ((*_matching)[u] != INVALID) {
deba@948
   193
          Node v = _graph.target((*_matching)[u]);
deba@948
   194
          _indeg->set(v, -1);
deba@948
   195
          Arc a = _graph.oppositeArc((*_matching)[u]);
deba@948
   196
          u = _graph.target((*_matching)[v]);
deba@948
   197
          _indeg->set(u, -1);
deba@948
   198
          _matching->set(v, a);
deba@948
   199
        }
deba@948
   200
      }
deba@948
   201
deba@948
   202
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
   203
        if ((*_indeg)[n] != 1) continue;
deba@948
   204
        _indeg->set(n, -1);
deba@948
   205
deba@948
   206
        int num = 1;
deba@948
   207
        Node u = _graph.target((*_matching)[n]);
deba@948
   208
        while (u != n) {
deba@948
   209
          _indeg->set(u, -1);
deba@948
   210
          u = _graph.target((*_matching)[u]);
deba@948
   211
          ++num;
deba@948
   212
        }
deba@948
   213
        if (num % 2 == 0 && num > 2) {
deba@948
   214
          Arc prev = _graph.oppositeArc((*_matching)[n]);
deba@948
   215
          Node v = _graph.target((*_matching)[n]);
deba@948
   216
          u = _graph.target((*_matching)[v]);
deba@948
   217
          _matching->set(v, prev);
deba@948
   218
          while (u != n) {
deba@948
   219
            prev = _graph.oppositeArc((*_matching)[u]);
deba@948
   220
            v = _graph.target((*_matching)[u]);
deba@948
   221
            u = _graph.target((*_matching)[v]);
deba@948
   222
            _matching->set(v, prev);
deba@948
   223
          }
deba@948
   224
        }
deba@948
   225
      }
deba@948
   226
    }
deba@948
   227
deba@948
   228
  public:
deba@948
   229
deba@948
   230
    typedef MaxFractionalMatching Create;
deba@948
   231
deba@948
   232
    ///\name Named Template Parameters
deba@948
   233
deba@948
   234
    ///@{
deba@948
   235
deba@948
   236
    template <typename T>
deba@948
   237
    struct SetMatchingMapTraits : public Traits {
deba@948
   238
      typedef T MatchingMap;
deba@948
   239
      static MatchingMap *createMatchingMap(const Graph&) {
deba@948
   240
        LEMON_ASSERT(false, "MatchingMap is not initialized");
deba@948
   241
        return 0; // ignore warnings
deba@948
   242
      }
deba@948
   243
    };
deba@948
   244
deba@948
   245
    /// \brief \ref named-templ-param "Named parameter" for setting
deba@948
   246
    /// MatchingMap type
deba@948
   247
    ///
deba@948
   248
    /// \ref named-templ-param "Named parameter" for setting MatchingMap
deba@948
   249
    /// type.
deba@948
   250
    template <typename T>
deba@948
   251
    struct SetMatchingMap
deba@948
   252
      : public MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > {
deba@948
   253
      typedef MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > Create;
deba@948
   254
    };
deba@948
   255
deba@948
   256
    template <typename T>
deba@948
   257
    struct SetElevatorTraits : public Traits {
deba@948
   258
      typedef T Elevator;
deba@948
   259
      static Elevator *createElevator(const Graph&, int) {
deba@948
   260
        LEMON_ASSERT(false, "Elevator is not initialized");
deba@948
   261
        return 0; // ignore warnings
deba@948
   262
      }
deba@948
   263
    };
deba@948
   264
deba@948
   265
    /// \brief \ref named-templ-param "Named parameter" for setting
deba@948
   266
    /// Elevator type
deba@948
   267
    ///
deba@948
   268
    /// \ref named-templ-param "Named parameter" for setting Elevator
deba@948
   269
    /// type. If this named parameter is used, then an external
deba@948
   270
    /// elevator object must be passed to the algorithm using the
deba@948
   271
    /// \ref elevator(Elevator&) "elevator()" function before calling
deba@948
   272
    /// \ref run() or \ref init().
deba@948
   273
    /// \sa SetStandardElevator
deba@948
   274
    template <typename T>
deba@948
   275
    struct SetElevator
deba@948
   276
      : public MaxFractionalMatching<Graph, SetElevatorTraits<T> > {
deba@948
   277
      typedef MaxFractionalMatching<Graph, SetElevatorTraits<T> > Create;
deba@948
   278
    };
deba@948
   279
deba@948
   280
    template <typename T>
deba@948
   281
    struct SetStandardElevatorTraits : public Traits {
deba@948
   282
      typedef T Elevator;
deba@948
   283
      static Elevator *createElevator(const Graph& graph, int max_level) {
deba@948
   284
        return new Elevator(graph, max_level);
deba@948
   285
      }
deba@948
   286
    };
deba@948
   287
deba@948
   288
    /// \brief \ref named-templ-param "Named parameter" for setting
deba@948
   289
    /// Elevator type with automatic allocation
deba@948
   290
    ///
deba@948
   291
    /// \ref named-templ-param "Named parameter" for setting Elevator
deba@948
   292
    /// type with automatic allocation.
deba@948
   293
    /// The Elevator should have standard constructor interface to be
deba@948
   294
    /// able to automatically created by the algorithm (i.e. the
deba@948
   295
    /// graph and the maximum level should be passed to it).
deba@948
   296
    /// However an external elevator object could also be passed to the
deba@948
   297
    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
deba@948
   298
    /// before calling \ref run() or \ref init().
deba@948
   299
    /// \sa SetElevator
deba@948
   300
    template <typename T>
deba@948
   301
    struct SetStandardElevator
deba@948
   302
      : public MaxFractionalMatching<Graph, SetStandardElevatorTraits<T> > {
deba@948
   303
      typedef MaxFractionalMatching<Graph,
deba@948
   304
                                    SetStandardElevatorTraits<T> > Create;
deba@948
   305
    };
deba@948
   306
deba@948
   307
    /// @}
deba@948
   308
deba@948
   309
  protected:
deba@948
   310
deba@948
   311
    MaxFractionalMatching() {}
deba@948
   312
deba@948
   313
  public:
deba@948
   314
deba@948
   315
    /// \brief Constructor
deba@948
   316
    ///
deba@948
   317
    /// Constructor.
deba@948
   318
    ///
deba@948
   319
    MaxFractionalMatching(const Graph &graph, bool allow_loops = true)
deba@948
   320
      : _graph(graph), _allow_loops(allow_loops),
deba@948
   321
        _local_matching(false), _matching(0),
deba@948
   322
        _local_level(false), _level(0),  _indeg(0)
deba@948
   323
    {}
deba@948
   324
deba@948
   325
    ~MaxFractionalMatching() {
deba@948
   326
      destroyStructures();
deba@948
   327
    }
deba@948
   328
deba@948
   329
    /// \brief Sets the matching map.
deba@948
   330
    ///
deba@948
   331
    /// Sets the matching map.
deba@948
   332
    /// If you don't use this function before calling \ref run() or
deba@948
   333
    /// \ref init(), an instance will be allocated automatically.
deba@948
   334
    /// The destructor deallocates this automatically allocated map,
deba@948
   335
    /// of course.
deba@948
   336
    /// \return <tt>(*this)</tt>
deba@948
   337
    MaxFractionalMatching& matchingMap(MatchingMap& map) {
deba@948
   338
      if (_local_matching) {
deba@948
   339
        delete _matching;
deba@948
   340
        _local_matching = false;
deba@948
   341
      }
deba@948
   342
      _matching = &map;
deba@948
   343
      return *this;
deba@948
   344
    }
deba@948
   345
deba@948
   346
    /// \brief Sets the elevator used by algorithm.
deba@948
   347
    ///
deba@948
   348
    /// Sets the elevator used by algorithm.
deba@948
   349
    /// If you don't use this function before calling \ref run() or
deba@948
   350
    /// \ref init(), an instance will be allocated automatically.
deba@948
   351
    /// The destructor deallocates this automatically allocated elevator,
deba@948
   352
    /// of course.
deba@948
   353
    /// \return <tt>(*this)</tt>
deba@948
   354
    MaxFractionalMatching& elevator(Elevator& elevator) {
deba@948
   355
      if (_local_level) {
deba@948
   356
        delete _level;
deba@948
   357
        _local_level = false;
deba@948
   358
      }
deba@948
   359
      _level = &elevator;
deba@948
   360
      return *this;
deba@948
   361
    }
deba@948
   362
deba@948
   363
    /// \brief Returns a const reference to the elevator.
deba@948
   364
    ///
deba@948
   365
    /// Returns a const reference to the elevator.
deba@948
   366
    ///
deba@948
   367
    /// \pre Either \ref run() or \ref init() must be called before
deba@948
   368
    /// using this function.
deba@948
   369
    const Elevator& elevator() const {
deba@948
   370
      return *_level;
deba@948
   371
    }
deba@948
   372
deba@948
   373
    /// \name Execution control
deba@948
   374
    /// The simplest way to execute the algorithm is to use one of the
deba@948
   375
    /// member functions called \c run(). \n
deba@948
   376
    /// If you need more control on the execution, first
deba@948
   377
    /// you must call \ref init() and then one variant of the start()
deba@948
   378
    /// member.
deba@948
   379
deba@948
   380
    /// @{
deba@948
   381
deba@948
   382
    /// \brief Initializes the internal data structures.
deba@948
   383
    ///
deba@948
   384
    /// Initializes the internal data structures and sets the initial
deba@948
   385
    /// matching.
deba@948
   386
    void init() {
deba@948
   387
      createStructures();
deba@948
   388
deba@948
   389
      _level->initStart();
deba@948
   390
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
   391
        _indeg->set(n, 0);
deba@948
   392
        _matching->set(n, INVALID);
deba@948
   393
        _level->initAddItem(n);
deba@948
   394
      }
deba@948
   395
      _level->initFinish();
deba@948
   396
deba@948
   397
      _empty_level = _node_num;
deba@948
   398
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
   399
        for (OutArcIt a(_graph, n); a != INVALID; ++a) {
deba@948
   400
          if (_graph.target(a) == n && !_allow_loops) continue;
deba@948
   401
          _matching->set(n, a);
deba@948
   402
          Node v = _graph.target((*_matching)[n]);
deba@948
   403
          _indeg->set(v, (*_indeg)[v] + 1);
deba@948
   404
          break;
deba@948
   405
        }
deba@948
   406
      }
deba@948
   407
deba@948
   408
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
   409
        if ((*_indeg)[n] == 0) {
deba@948
   410
          _level->activate(n);
deba@948
   411
        }
deba@948
   412
      }
deba@948
   413
    }
deba@948
   414
deba@948
   415
    /// \brief Starts the algorithm and computes a fractional matching
deba@948
   416
    ///
deba@948
   417
    /// The algorithm computes a maximum fractional matching.
deba@948
   418
    ///
deba@948
   419
    /// \param postprocess The algorithm computes first a matching
deba@948
   420
    /// which is a union of a matching with one value edges, cycles
deba@948
   421
    /// with half value edges and even length paths with half value
deba@948
   422
    /// edges. If the parameter is true, then after the push-relabel
deba@948
   423
    /// algorithm it postprocesses the matching to contain only
deba@948
   424
    /// matching edges and half value odd cycles.
deba@948
   425
    void start(bool postprocess = true) {
deba@948
   426
      Node n;
deba@948
   427
      while ((n = _level->highestActive()) != INVALID) {
deba@948
   428
        int level = _level->highestActiveLevel();
deba@948
   429
        int new_level = _level->maxLevel();
deba@948
   430
        for (InArcIt a(_graph, n); a != INVALID; ++a) {
deba@948
   431
          Node u = _graph.source(a);
deba@948
   432
          if (n == u && !_allow_loops) continue;
deba@948
   433
          Node v = _graph.target((*_matching)[u]);
deba@948
   434
          if ((*_level)[v] < level) {
deba@948
   435
            _indeg->set(v, (*_indeg)[v] - 1);
deba@948
   436
            if ((*_indeg)[v] == 0) {
deba@948
   437
              _level->activate(v);
deba@948
   438
            }
deba@948
   439
            _matching->set(u, a);
deba@948
   440
            _indeg->set(n, (*_indeg)[n] + 1);
deba@948
   441
            _level->deactivate(n);
deba@948
   442
            goto no_more_push;
deba@948
   443
          } else if (new_level > (*_level)[v]) {
deba@948
   444
            new_level = (*_level)[v];
deba@948
   445
          }
deba@948
   446
        }
deba@948
   447
deba@948
   448
        if (new_level + 1 < _level->maxLevel()) {
deba@948
   449
          _level->liftHighestActive(new_level + 1);
deba@948
   450
        } else {
deba@948
   451
          _level->liftHighestActiveToTop();
deba@948
   452
        }
deba@948
   453
        if (_level->emptyLevel(level)) {
deba@948
   454
          _level->liftToTop(level);
deba@948
   455
        }
deba@948
   456
      no_more_push:
deba@948
   457
        ;
deba@948
   458
      }
deba@948
   459
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
   460
        if ((*_matching)[n] == INVALID) continue;
deba@948
   461
        Node u = _graph.target((*_matching)[n]);
deba@948
   462
        if ((*_indeg)[u] > 1) {
deba@948
   463
          _indeg->set(u, (*_indeg)[u] - 1);
deba@948
   464
          _matching->set(n, INVALID);
deba@948
   465
        }
deba@948
   466
      }
deba@948
   467
      if (postprocess) {
deba@948
   468
        postprocessing();
deba@948
   469
      }
deba@948
   470
    }
deba@948
   471
deba@948
   472
    /// \brief Starts the algorithm and computes a perfect fractional
deba@948
   473
    /// matching
deba@948
   474
    ///
deba@948
   475
    /// The algorithm computes a perfect fractional matching. If it
deba@948
   476
    /// does not exists, then the algorithm returns false and the
deba@948
   477
    /// matching is undefined and the barrier.
deba@948
   478
    ///
deba@948
   479
    /// \param postprocess The algorithm computes first a matching
deba@948
   480
    /// which is a union of a matching with one value edges, cycles
deba@948
   481
    /// with half value edges and even length paths with half value
deba@948
   482
    /// edges. If the parameter is true, then after the push-relabel
deba@948
   483
    /// algorithm it postprocesses the matching to contain only
deba@948
   484
    /// matching edges and half value odd cycles.
deba@948
   485
    bool startPerfect(bool postprocess = true) {
deba@948
   486
      Node n;
deba@948
   487
      while ((n = _level->highestActive()) != INVALID) {
deba@948
   488
        int level = _level->highestActiveLevel();
deba@948
   489
        int new_level = _level->maxLevel();
deba@948
   490
        for (InArcIt a(_graph, n); a != INVALID; ++a) {
deba@948
   491
          Node u = _graph.source(a);
deba@948
   492
          if (n == u && !_allow_loops) continue;
deba@948
   493
          Node v = _graph.target((*_matching)[u]);
deba@948
   494
          if ((*_level)[v] < level) {
deba@948
   495
            _indeg->set(v, (*_indeg)[v] - 1);
deba@948
   496
            if ((*_indeg)[v] == 0) {
deba@948
   497
              _level->activate(v);
deba@948
   498
            }
deba@948
   499
            _matching->set(u, a);
deba@948
   500
            _indeg->set(n, (*_indeg)[n] + 1);
deba@948
   501
            _level->deactivate(n);
deba@948
   502
            goto no_more_push;
deba@948
   503
          } else if (new_level > (*_level)[v]) {
deba@948
   504
            new_level = (*_level)[v];
deba@948
   505
          }
deba@948
   506
        }
deba@948
   507
deba@948
   508
        if (new_level + 1 < _level->maxLevel()) {
deba@948
   509
          _level->liftHighestActive(new_level + 1);
deba@948
   510
        } else {
deba@948
   511
          _level->liftHighestActiveToTop();
deba@948
   512
          _empty_level = _level->maxLevel() - 1;
deba@948
   513
          return false;
deba@948
   514
        }
deba@948
   515
        if (_level->emptyLevel(level)) {
deba@948
   516
          _level->liftToTop(level);
deba@948
   517
          _empty_level = level;
deba@948
   518
          return false;
deba@948
   519
        }
deba@948
   520
      no_more_push:
deba@948
   521
        ;
deba@948
   522
      }
deba@948
   523
      if (postprocess) {
deba@948
   524
        postprocessing();
deba@948
   525
      }
deba@948
   526
      return true;
deba@948
   527
    }
deba@948
   528
deba@948
   529
    /// \brief Runs the algorithm
deba@948
   530
    ///
deba@948
   531
    /// Just a shortcut for the next code:
deba@948
   532
    ///\code
deba@948
   533
    /// init();
deba@948
   534
    /// start();
deba@948
   535
    ///\endcode
deba@948
   536
    void run(bool postprocess = true) {
deba@948
   537
      init();
deba@948
   538
      start(postprocess);
deba@948
   539
    }
deba@948
   540
deba@948
   541
    /// \brief Runs the algorithm to find a perfect fractional matching
deba@948
   542
    ///
deba@948
   543
    /// Just a shortcut for the next code:
deba@948
   544
    ///\code
deba@948
   545
    /// init();
deba@948
   546
    /// startPerfect();
deba@948
   547
    ///\endcode
deba@948
   548
    bool runPerfect(bool postprocess = true) {
deba@948
   549
      init();
deba@948
   550
      return startPerfect(postprocess);
deba@948
   551
    }
deba@948
   552
deba@948
   553
    ///@}
deba@948
   554
deba@948
   555
    /// \name Query Functions
deba@948
   556
    /// The result of the %Matching algorithm can be obtained using these
deba@948
   557
    /// functions.\n
deba@948
   558
    /// Before the use of these functions,
deba@948
   559
    /// either run() or start() must be called.
deba@948
   560
    ///@{
deba@948
   561
deba@948
   562
deba@948
   563
    /// \brief Return the number of covered nodes in the matching.
deba@948
   564
    ///
deba@948
   565
    /// This function returns the number of covered nodes in the matching.
deba@948
   566
    ///
deba@948
   567
    /// \pre Either run() or start() must be called before using this function.
deba@948
   568
    int matchingSize() const {
deba@948
   569
      int num = 0;
deba@948
   570
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
   571
        if ((*_matching)[n] != INVALID) {
deba@948
   572
          ++num;
deba@948
   573
        }
deba@948
   574
      }
deba@948
   575
      return num;
deba@948
   576
    }
deba@948
   577
deba@948
   578
    /// \brief Returns a const reference to the matching map.
deba@948
   579
    ///
deba@948
   580
    /// Returns a const reference to the node map storing the found
deba@948
   581
    /// fractional matching. This method can be called after
deba@948
   582
    /// running the algorithm.
deba@948
   583
    ///
deba@948
   584
    /// \pre Either \ref run() or \ref init() must be called before
deba@948
   585
    /// using this function.
deba@948
   586
    const MatchingMap& matchingMap() const {
deba@948
   587
      return *_matching;
deba@948
   588
    }
deba@948
   589
deba@948
   590
    /// \brief Return \c true if the given edge is in the matching.
deba@948
   591
    ///
deba@948
   592
    /// This function returns \c true if the given edge is in the
deba@948
   593
    /// found matching. The result is scaled by \ref primalScale
deba@948
   594
    /// "primal scale".
deba@948
   595
    ///
deba@948
   596
    /// \pre Either run() or start() must be called before using this function.
deba@948
   597
    int matching(const Edge& edge) const {
deba@948
   598
      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0) +
deba@948
   599
        (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
deba@948
   600
    }
deba@948
   601
deba@948
   602
    /// \brief Return the fractional matching arc (or edge) incident
deba@948
   603
    /// to the given node.
deba@948
   604
    ///
deba@948
   605
    /// This function returns one of the fractional matching arc (or
deba@948
   606
    /// edge) incident to the given node in the found matching or \c
deba@948
   607
    /// INVALID if the node is not covered by the matching or if the
deba@948
   608
    /// node is on an odd length cycle then it is the successor edge
deba@948
   609
    /// on the cycle.
deba@948
   610
    ///
deba@948
   611
    /// \pre Either run() or start() must be called before using this function.
deba@948
   612
    Arc matching(const Node& node) const {
deba@948
   613
      return (*_matching)[node];
deba@948
   614
    }
deba@948
   615
deba@948
   616
    /// \brief Returns true if the node is in the barrier
deba@948
   617
    ///
deba@948
   618
    /// The barrier is a subset of the nodes. If the nodes in the
deba@948
   619
    /// barrier have less adjacent nodes than the size of the barrier,
deba@948
   620
    /// then at least as much nodes cannot be covered as the
deba@948
   621
    /// difference of the two subsets.
deba@948
   622
    bool barrier(const Node& node) const {
deba@948
   623
      return (*_level)[node] >= _empty_level;
deba@948
   624
    }
deba@948
   625
deba@948
   626
    /// @}
deba@948
   627
deba@948
   628
  };
deba@948
   629
deba@948
   630
  /// \ingroup matching
deba@948
   631
  ///
deba@948
   632
  /// \brief Weighted fractional matching in general graphs
deba@948
   633
  ///
deba@948
   634
  /// This class provides an efficient implementation of fractional
deba@950
   635
  /// matching algorithm. The implementation uses priority queues and
deba@950
   636
  /// provides \f$O(nm\log n)\f$ time complexity.
deba@948
   637
  ///
deba@948
   638
  /// The maximum weighted fractional matching is a relaxation of the
deba@948
   639
  /// maximum weighted matching problem where the odd set constraints
deba@948
   640
  /// are omitted.
deba@948
   641
  /// It can be formulated with the following linear program.
deba@948
   642
  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
deba@948
   643
  /// \f[x_e \ge 0\quad \forall e\in E\f]
deba@948
   644
  /// \f[\max \sum_{e\in E}x_ew_e\f]
deba@948
   645
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
deba@948
   646
  /// \f$X\f$. The result must be the union of a matching with one
deba@948
   647
  /// value edges and a set of odd length cycles with half value edges.
deba@948
   648
  ///
deba@948
   649
  /// The algorithm calculates an optimal fractional matching and a
deba@948
   650
  /// proof of the optimality. The solution of the dual problem can be
deba@948
   651
  /// used to check the result of the algorithm. The dual linear
deba@948
   652
  /// problem is the following.
deba@948
   653
  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
deba@948
   654
  /// \f[y_u \ge 0 \quad \forall u \in V\f]
deba@950
   655
  /// \f[\min \sum_{u \in V}y_u \f]
deba@948
   656
  ///
deba@948
   657
  /// The algorithm can be executed with the run() function.
deba@948
   658
  /// After it the matching (the primal solution) and the dual solution
deba@948
   659
  /// can be obtained using the query functions.
deba@948
   660
  ///
deba@948
   661
  /// If the value type is integer, then the primal and the dual
deba@948
   662
  /// solutions are multiplied by
deba@950
   663
  /// \ref MaxWeightedFractionalMatching::primalScale "2" and
deba@950
   664
  /// \ref MaxWeightedFractionalMatching::dualScale "4" respectively.
deba@948
   665
  ///
deba@948
   666
  /// \tparam GR The undirected graph type the algorithm runs on.
deba@948
   667
  /// \tparam WM The type edge weight map. The default type is
deba@948
   668
  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
deba@948
   669
#ifdef DOXYGEN
deba@948
   670
  template <typename GR, typename WM>
deba@948
   671
#else
deba@948
   672
  template <typename GR,
deba@948
   673
            typename WM = typename GR::template EdgeMap<int> >
deba@948
   674
#endif
deba@948
   675
  class MaxWeightedFractionalMatching {
deba@948
   676
  public:
deba@948
   677
deba@948
   678
    /// The graph type of the algorithm
deba@948
   679
    typedef GR Graph;
deba@948
   680
    /// The type of the edge weight map
deba@948
   681
    typedef WM WeightMap;
deba@948
   682
    /// The value type of the edge weights
deba@948
   683
    typedef typename WeightMap::Value Value;
deba@948
   684
deba@948
   685
    /// The type of the matching map
deba@948
   686
    typedef typename Graph::template NodeMap<typename Graph::Arc>
deba@948
   687
    MatchingMap;
deba@948
   688
deba@948
   689
    /// \brief Scaling factor for primal solution
deba@948
   690
    ///
deba@948
   691
    /// Scaling factor for primal solution. It is equal to 2 or 1
deba@948
   692
    /// according to the value type.
deba@948
   693
    static const int primalScale =
deba@948
   694
      std::numeric_limits<Value>::is_integer ? 2 : 1;
deba@948
   695
deba@948
   696
    /// \brief Scaling factor for dual solution
deba@948
   697
    ///
deba@948
   698
    /// Scaling factor for dual solution. It is equal to 4 or 1
deba@948
   699
    /// according to the value type.
deba@948
   700
    static const int dualScale =
deba@948
   701
      std::numeric_limits<Value>::is_integer ? 4 : 1;
deba@948
   702
deba@948
   703
  private:
deba@948
   704
deba@948
   705
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
deba@948
   706
deba@948
   707
    typedef typename Graph::template NodeMap<Value> NodePotential;
deba@948
   708
deba@948
   709
    const Graph& _graph;
deba@948
   710
    const WeightMap& _weight;
deba@948
   711
deba@948
   712
    MatchingMap* _matching;
deba@948
   713
    NodePotential* _node_potential;
deba@948
   714
deba@948
   715
    int _node_num;
deba@948
   716
    bool _allow_loops;
deba@948
   717
deba@948
   718
    enum Status {
deba@948
   719
      EVEN = -1, MATCHED = 0, ODD = 1
deba@948
   720
    };
deba@948
   721
deba@948
   722
    typedef typename Graph::template NodeMap<Status> StatusMap;
deba@948
   723
    StatusMap* _status;
deba@948
   724
deba@948
   725
    typedef typename Graph::template NodeMap<Arc> PredMap;
deba@948
   726
    PredMap* _pred;
deba@948
   727
deba@948
   728
    typedef ExtendFindEnum<IntNodeMap> TreeSet;
deba@948
   729
deba@948
   730
    IntNodeMap *_tree_set_index;
deba@948
   731
    TreeSet *_tree_set;
deba@948
   732
deba@948
   733
    IntNodeMap *_delta1_index;
deba@948
   734
    BinHeap<Value, IntNodeMap> *_delta1;
deba@948
   735
deba@948
   736
    IntNodeMap *_delta2_index;
deba@948
   737
    BinHeap<Value, IntNodeMap> *_delta2;
deba@948
   738
deba@948
   739
    IntEdgeMap *_delta3_index;
deba@948
   740
    BinHeap<Value, IntEdgeMap> *_delta3;
deba@948
   741
deba@948
   742
    Value _delta_sum;
deba@948
   743
deba@948
   744
    void createStructures() {
deba@948
   745
      _node_num = countNodes(_graph);
deba@948
   746
deba@948
   747
      if (!_matching) {
deba@948
   748
        _matching = new MatchingMap(_graph);
deba@948
   749
      }
deba@948
   750
      if (!_node_potential) {
deba@948
   751
        _node_potential = new NodePotential(_graph);
deba@948
   752
      }
deba@948
   753
      if (!_status) {
deba@948
   754
        _status = new StatusMap(_graph);
deba@948
   755
      }
deba@948
   756
      if (!_pred) {
deba@948
   757
        _pred = new PredMap(_graph);
deba@948
   758
      }
deba@948
   759
      if (!_tree_set) {
deba@948
   760
        _tree_set_index = new IntNodeMap(_graph);
deba@948
   761
        _tree_set = new TreeSet(*_tree_set_index);
deba@948
   762
      }
deba@948
   763
      if (!_delta1) {
deba@948
   764
        _delta1_index = new IntNodeMap(_graph);
deba@948
   765
        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
deba@948
   766
      }
deba@948
   767
      if (!_delta2) {
deba@948
   768
        _delta2_index = new IntNodeMap(_graph);
deba@948
   769
        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
deba@948
   770
      }
deba@948
   771
      if (!_delta3) {
deba@948
   772
        _delta3_index = new IntEdgeMap(_graph);
deba@948
   773
        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
deba@948
   774
      }
deba@948
   775
    }
deba@948
   776
deba@948
   777
    void destroyStructures() {
deba@948
   778
      if (_matching) {
deba@948
   779
        delete _matching;
deba@948
   780
      }
deba@948
   781
      if (_node_potential) {
deba@948
   782
        delete _node_potential;
deba@948
   783
      }
deba@948
   784
      if (_status) {
deba@948
   785
        delete _status;
deba@948
   786
      }
deba@948
   787
      if (_pred) {
deba@948
   788
        delete _pred;
deba@948
   789
      }
deba@948
   790
      if (_tree_set) {
deba@948
   791
        delete _tree_set_index;
deba@948
   792
        delete _tree_set;
deba@948
   793
      }
deba@948
   794
      if (_delta1) {
deba@948
   795
        delete _delta1_index;
deba@948
   796
        delete _delta1;
deba@948
   797
      }
deba@948
   798
      if (_delta2) {
deba@948
   799
        delete _delta2_index;
deba@948
   800
        delete _delta2;
deba@948
   801
      }
deba@948
   802
      if (_delta3) {
deba@948
   803
        delete _delta3_index;
deba@948
   804
        delete _delta3;
deba@948
   805
      }
deba@948
   806
    }
deba@948
   807
deba@948
   808
    void matchedToEven(Node node, int tree) {
deba@948
   809
      _tree_set->insert(node, tree);
deba@948
   810
      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
deba@948
   811
      _delta1->push(node, (*_node_potential)[node]);
deba@948
   812
deba@948
   813
      if (_delta2->state(node) == _delta2->IN_HEAP) {
deba@948
   814
        _delta2->erase(node);
deba@948
   815
      }
deba@948
   816
deba@948
   817
      for (InArcIt a(_graph, node); a != INVALID; ++a) {
deba@948
   818
        Node v = _graph.source(a);
deba@948
   819
        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
deba@948
   820
          dualScale * _weight[a];
deba@948
   821
        if (node == v) {
deba@948
   822
          if (_allow_loops && _graph.direction(a)) {
deba@948
   823
            _delta3->push(a, rw / 2);
deba@948
   824
          }
deba@948
   825
        } else if ((*_status)[v] == EVEN) {
deba@948
   826
          _delta3->push(a, rw / 2);
deba@948
   827
        } else if ((*_status)[v] == MATCHED) {
deba@948
   828
          if (_delta2->state(v) != _delta2->IN_HEAP) {
deba@948
   829
            _pred->set(v, a);
deba@948
   830
            _delta2->push(v, rw);
deba@948
   831
          } else if ((*_delta2)[v] > rw) {
deba@948
   832
            _pred->set(v, a);
deba@948
   833
            _delta2->decrease(v, rw);
deba@948
   834
          }
deba@948
   835
        }
deba@948
   836
      }
deba@948
   837
    }
deba@948
   838
deba@948
   839
    void matchedToOdd(Node node, int tree) {
deba@948
   840
      _tree_set->insert(node, tree);
deba@948
   841
      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
deba@948
   842
deba@948
   843
      if (_delta2->state(node) == _delta2->IN_HEAP) {
deba@948
   844
        _delta2->erase(node);
deba@948
   845
      }
deba@948
   846
    }
deba@948
   847
deba@948
   848
    void evenToMatched(Node node, int tree) {
deba@948
   849
      _delta1->erase(node);
deba@948
   850
      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
deba@948
   851
      Arc min = INVALID;
deba@948
   852
      Value minrw = std::numeric_limits<Value>::max();
deba@948
   853
      for (InArcIt a(_graph, node); a != INVALID; ++a) {
deba@948
   854
        Node v = _graph.source(a);
deba@948
   855
        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
deba@948
   856
          dualScale * _weight[a];
deba@948
   857
deba@948
   858
        if (node == v) {
deba@948
   859
          if (_allow_loops && _graph.direction(a)) {
deba@948
   860
            _delta3->erase(a);
deba@948
   861
          }
deba@948
   862
        } else if ((*_status)[v] == EVEN) {
deba@948
   863
          _delta3->erase(a);
deba@948
   864
          if (minrw > rw) {
deba@948
   865
            min = _graph.oppositeArc(a);
deba@948
   866
            minrw = rw;
deba@948
   867
          }
deba@948
   868
        } else if ((*_status)[v]  == MATCHED) {
deba@948
   869
          if ((*_pred)[v] == a) {
deba@948
   870
            Arc mina = INVALID;
deba@948
   871
            Value minrwa = std::numeric_limits<Value>::max();
deba@948
   872
            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
deba@948
   873
              Node va = _graph.target(aa);
deba@948
   874
              if ((*_status)[va] != EVEN ||
deba@948
   875
                  _tree_set->find(va) == tree) continue;
deba@948
   876
              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
deba@948
   877
                dualScale * _weight[aa];
deba@948
   878
              if (minrwa > rwa) {
deba@948
   879
                minrwa = rwa;
deba@948
   880
                mina = aa;
deba@948
   881
              }
deba@948
   882
            }
deba@948
   883
            if (mina != INVALID) {
deba@948
   884
              _pred->set(v, mina);
deba@948
   885
              _delta2->increase(v, minrwa);
deba@948
   886
            } else {
deba@948
   887
              _pred->set(v, INVALID);
deba@948
   888
              _delta2->erase(v);
deba@948
   889
            }
deba@948
   890
          }
deba@948
   891
        }
deba@948
   892
      }
deba@948
   893
      if (min != INVALID) {
deba@948
   894
        _pred->set(node, min);
deba@948
   895
        _delta2->push(node, minrw);
deba@948
   896
      } else {
deba@948
   897
        _pred->set(node, INVALID);
deba@948
   898
      }
deba@948
   899
    }
deba@948
   900
deba@948
   901
    void oddToMatched(Node node) {
deba@948
   902
      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
deba@948
   903
      Arc min = INVALID;
deba@948
   904
      Value minrw = std::numeric_limits<Value>::max();
deba@948
   905
      for (InArcIt a(_graph, node); a != INVALID; ++a) {
deba@948
   906
        Node v = _graph.source(a);
deba@948
   907
        if ((*_status)[v] != EVEN) continue;
deba@948
   908
        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
deba@948
   909
          dualScale * _weight[a];
deba@948
   910
deba@948
   911
        if (minrw > rw) {
deba@948
   912
          min = _graph.oppositeArc(a);
deba@948
   913
          minrw = rw;
deba@948
   914
        }
deba@948
   915
      }
deba@948
   916
      if (min != INVALID) {
deba@948
   917
        _pred->set(node, min);
deba@948
   918
        _delta2->push(node, minrw);
deba@948
   919
      } else {
deba@948
   920
        _pred->set(node, INVALID);
deba@948
   921
      }
deba@948
   922
    }
deba@948
   923
deba@948
   924
    void alternatePath(Node even, int tree) {
deba@948
   925
      Node odd;
deba@948
   926
deba@948
   927
      _status->set(even, MATCHED);
deba@948
   928
      evenToMatched(even, tree);
deba@948
   929
deba@948
   930
      Arc prev = (*_matching)[even];
deba@948
   931
      while (prev != INVALID) {
deba@948
   932
        odd = _graph.target(prev);
deba@948
   933
        even = _graph.target((*_pred)[odd]);
deba@948
   934
        _matching->set(odd, (*_pred)[odd]);
deba@948
   935
        _status->set(odd, MATCHED);
deba@948
   936
        oddToMatched(odd);
deba@948
   937
deba@948
   938
        prev = (*_matching)[even];
deba@948
   939
        _status->set(even, MATCHED);
deba@948
   940
        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
deba@948
   941
        evenToMatched(even, tree);
deba@948
   942
      }
deba@948
   943
    }
deba@948
   944
deba@948
   945
    void destroyTree(int tree) {
deba@948
   946
      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
deba@948
   947
        if ((*_status)[n] == EVEN) {
deba@948
   948
          _status->set(n, MATCHED);
deba@948
   949
          evenToMatched(n, tree);
deba@948
   950
        } else if ((*_status)[n] == ODD) {
deba@948
   951
          _status->set(n, MATCHED);
deba@948
   952
          oddToMatched(n);
deba@948
   953
        }
deba@948
   954
      }
deba@948
   955
      _tree_set->eraseClass(tree);
deba@948
   956
    }
deba@948
   957
deba@948
   958
deba@948
   959
    void unmatchNode(const Node& node) {
deba@948
   960
      int tree = _tree_set->find(node);
deba@948
   961
deba@948
   962
      alternatePath(node, tree);
deba@948
   963
      destroyTree(tree);
deba@948
   964
deba@948
   965
      _matching->set(node, INVALID);
deba@948
   966
    }
deba@948
   967
deba@948
   968
deba@948
   969
    void augmentOnEdge(const Edge& edge) {
deba@948
   970
      Node left = _graph.u(edge);
deba@948
   971
      int left_tree = _tree_set->find(left);
deba@948
   972
deba@948
   973
      alternatePath(left, left_tree);
deba@948
   974
      destroyTree(left_tree);
deba@948
   975
      _matching->set(left, _graph.direct(edge, true));
deba@948
   976
deba@948
   977
      Node right = _graph.v(edge);
deba@948
   978
      int right_tree = _tree_set->find(right);
deba@948
   979
deba@948
   980
      alternatePath(right, right_tree);
deba@948
   981
      destroyTree(right_tree);
deba@948
   982
      _matching->set(right, _graph.direct(edge, false));
deba@948
   983
    }
deba@948
   984
deba@948
   985
    void augmentOnArc(const Arc& arc) {
deba@948
   986
      Node left = _graph.source(arc);
deba@948
   987
      _status->set(left, MATCHED);
deba@948
   988
      _matching->set(left, arc);
deba@948
   989
      _pred->set(left, arc);
deba@948
   990
deba@948
   991
      Node right = _graph.target(arc);
deba@948
   992
      int right_tree = _tree_set->find(right);
deba@948
   993
deba@948
   994
      alternatePath(right, right_tree);
deba@948
   995
      destroyTree(right_tree);
deba@948
   996
      _matching->set(right, _graph.oppositeArc(arc));
deba@948
   997
    }
deba@948
   998
deba@948
   999
    void extendOnArc(const Arc& arc) {
deba@948
  1000
      Node base = _graph.target(arc);
deba@948
  1001
      int tree = _tree_set->find(base);
deba@948
  1002
deba@948
  1003
      Node odd = _graph.source(arc);
deba@948
  1004
      _tree_set->insert(odd, tree);
deba@948
  1005
      _status->set(odd, ODD);
deba@948
  1006
      matchedToOdd(odd, tree);
deba@948
  1007
      _pred->set(odd, arc);
deba@948
  1008
deba@948
  1009
      Node even = _graph.target((*_matching)[odd]);
deba@948
  1010
      _tree_set->insert(even, tree);
deba@948
  1011
      _status->set(even, EVEN);
deba@948
  1012
      matchedToEven(even, tree);
deba@948
  1013
    }
deba@948
  1014
deba@948
  1015
    void cycleOnEdge(const Edge& edge, int tree) {
deba@948
  1016
      Node nca = INVALID;
deba@948
  1017
      std::vector<Node> left_path, right_path;
deba@948
  1018
deba@948
  1019
      {
deba@948
  1020
        std::set<Node> left_set, right_set;
deba@948
  1021
        Node left = _graph.u(edge);
deba@948
  1022
        left_path.push_back(left);
deba@948
  1023
        left_set.insert(left);
deba@948
  1024
deba@948
  1025
        Node right = _graph.v(edge);
deba@948
  1026
        right_path.push_back(right);
deba@948
  1027
        right_set.insert(right);
deba@948
  1028
deba@948
  1029
        while (true) {
deba@948
  1030
deba@948
  1031
          if (left_set.find(right) != left_set.end()) {
deba@948
  1032
            nca = right;
deba@948
  1033
            break;
deba@948
  1034
          }
deba@948
  1035
deba@948
  1036
          if ((*_matching)[left] == INVALID) break;
deba@948
  1037
deba@948
  1038
          left = _graph.target((*_matching)[left]);
deba@948
  1039
          left_path.push_back(left);
deba@948
  1040
          left = _graph.target((*_pred)[left]);
deba@948
  1041
          left_path.push_back(left);
deba@948
  1042
deba@948
  1043
          left_set.insert(left);
deba@948
  1044
deba@948
  1045
          if (right_set.find(left) != right_set.end()) {
deba@948
  1046
            nca = left;
deba@948
  1047
            break;
deba@948
  1048
          }
deba@948
  1049
deba@948
  1050
          if ((*_matching)[right] == INVALID) break;
deba@948
  1051
deba@948
  1052
          right = _graph.target((*_matching)[right]);
deba@948
  1053
          right_path.push_back(right);
deba@948
  1054
          right = _graph.target((*_pred)[right]);
deba@948
  1055
          right_path.push_back(right);
deba@948
  1056
deba@948
  1057
          right_set.insert(right);
deba@948
  1058
deba@948
  1059
        }
deba@948
  1060
deba@948
  1061
        if (nca == INVALID) {
deba@948
  1062
          if ((*_matching)[left] == INVALID) {
deba@948
  1063
            nca = right;
deba@948
  1064
            while (left_set.find(nca) == left_set.end()) {
deba@948
  1065
              nca = _graph.target((*_matching)[nca]);
deba@948
  1066
              right_path.push_back(nca);
deba@948
  1067
              nca = _graph.target((*_pred)[nca]);
deba@948
  1068
              right_path.push_back(nca);
deba@948
  1069
            }
deba@948
  1070
          } else {
deba@948
  1071
            nca = left;
deba@948
  1072
            while (right_set.find(nca) == right_set.end()) {
deba@948
  1073
              nca = _graph.target((*_matching)[nca]);
deba@948
  1074
              left_path.push_back(nca);
deba@948
  1075
              nca = _graph.target((*_pred)[nca]);
deba@948
  1076
              left_path.push_back(nca);
deba@948
  1077
            }
deba@948
  1078
          }
deba@948
  1079
        }
deba@948
  1080
      }
deba@948
  1081
deba@948
  1082
      alternatePath(nca, tree);
deba@948
  1083
      Arc prev;
deba@948
  1084
deba@948
  1085
      prev = _graph.direct(edge, true);
deba@948
  1086
      for (int i = 0; left_path[i] != nca; i += 2) {
deba@948
  1087
        _matching->set(left_path[i], prev);
deba@948
  1088
        _status->set(left_path[i], MATCHED);
deba@948
  1089
        evenToMatched(left_path[i], tree);
deba@948
  1090
deba@948
  1091
        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
deba@948
  1092
        _status->set(left_path[i + 1], MATCHED);
deba@948
  1093
        oddToMatched(left_path[i + 1]);
deba@948
  1094
      }
deba@948
  1095
      _matching->set(nca, prev);
deba@948
  1096
deba@948
  1097
      for (int i = 0; right_path[i] != nca; i += 2) {
deba@948
  1098
        _status->set(right_path[i], MATCHED);
deba@948
  1099
        evenToMatched(right_path[i], tree);
deba@948
  1100
deba@948
  1101
        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
deba@948
  1102
        _status->set(right_path[i + 1], MATCHED);
deba@948
  1103
        oddToMatched(right_path[i + 1]);
deba@948
  1104
      }
deba@948
  1105
deba@948
  1106
      destroyTree(tree);
deba@948
  1107
    }
deba@948
  1108
deba@948
  1109
    void extractCycle(const Arc &arc) {
deba@948
  1110
      Node left = _graph.source(arc);
deba@948
  1111
      Node odd = _graph.target((*_matching)[left]);
deba@948
  1112
      Arc prev;
deba@948
  1113
      while (odd != left) {
deba@948
  1114
        Node even = _graph.target((*_matching)[odd]);
deba@948
  1115
        prev = (*_matching)[odd];
deba@948
  1116
        odd = _graph.target((*_matching)[even]);
deba@948
  1117
        _matching->set(even, _graph.oppositeArc(prev));
deba@948
  1118
      }
deba@948
  1119
      _matching->set(left, arc);
deba@948
  1120
deba@948
  1121
      Node right = _graph.target(arc);
deba@948
  1122
      int right_tree = _tree_set->find(right);
deba@948
  1123
      alternatePath(right, right_tree);
deba@948
  1124
      destroyTree(right_tree);
deba@948
  1125
      _matching->set(right, _graph.oppositeArc(arc));
deba@948
  1126
    }
deba@948
  1127
deba@948
  1128
  public:
deba@948
  1129
deba@948
  1130
    /// \brief Constructor
deba@948
  1131
    ///
deba@948
  1132
    /// Constructor.
deba@948
  1133
    MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
deba@948
  1134
                                  bool allow_loops = true)
deba@948
  1135
      : _graph(graph), _weight(weight), _matching(0),
deba@948
  1136
      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
deba@948
  1137
      _status(0),  _pred(0),
deba@948
  1138
      _tree_set_index(0), _tree_set(0),
deba@948
  1139
deba@948
  1140
      _delta1_index(0), _delta1(0),
deba@948
  1141
      _delta2_index(0), _delta2(0),
deba@948
  1142
      _delta3_index(0), _delta3(0),
deba@948
  1143
deba@948
  1144
      _delta_sum() {}
deba@948
  1145
deba@948
  1146
    ~MaxWeightedFractionalMatching() {
deba@948
  1147
      destroyStructures();
deba@948
  1148
    }
deba@948
  1149
deba@948
  1150
    /// \name Execution Control
deba@948
  1151
    /// The simplest way to execute the algorithm is to use the
deba@948
  1152
    /// \ref run() member function.
deba@948
  1153
deba@948
  1154
    ///@{
deba@948
  1155
deba@948
  1156
    /// \brief Initialize the algorithm
deba@948
  1157
    ///
deba@948
  1158
    /// This function initializes the algorithm.
deba@948
  1159
    void init() {
deba@948
  1160
      createStructures();
deba@948
  1161
deba@948
  1162
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
  1163
        (*_delta1_index)[n] = _delta1->PRE_HEAP;
deba@948
  1164
        (*_delta2_index)[n] = _delta2->PRE_HEAP;
deba@948
  1165
      }
deba@948
  1166
      for (EdgeIt e(_graph); e != INVALID; ++e) {
deba@948
  1167
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
deba@948
  1168
      }
deba@948
  1169
deba@948
  1170
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
  1171
        Value max = 0;
deba@948
  1172
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
deba@948
  1173
          if (_graph.target(e) == n && !_allow_loops) continue;
deba@948
  1174
          if ((dualScale * _weight[e]) / 2 > max) {
deba@948
  1175
            max = (dualScale * _weight[e]) / 2;
deba@948
  1176
          }
deba@948
  1177
        }
deba@948
  1178
        _node_potential->set(n, max);
deba@948
  1179
        _delta1->push(n, max);
deba@948
  1180
deba@948
  1181
        _tree_set->insert(n);
deba@948
  1182
deba@948
  1183
        _matching->set(n, INVALID);
deba@948
  1184
        _status->set(n, EVEN);
deba@948
  1185
      }
deba@948
  1186
deba@948
  1187
      for (EdgeIt e(_graph); e != INVALID; ++e) {
deba@948
  1188
        Node left = _graph.u(e);
deba@948
  1189
        Node right = _graph.v(e);
deba@948
  1190
        if (left == right && !_allow_loops) continue;
deba@948
  1191
        _delta3->push(e, ((*_node_potential)[left] +
deba@948
  1192
                          (*_node_potential)[right] -
deba@948
  1193
                          dualScale * _weight[e]) / 2);
deba@948
  1194
      }
deba@948
  1195
    }
deba@948
  1196
deba@948
  1197
    /// \brief Start the algorithm
deba@948
  1198
    ///
deba@948
  1199
    /// This function starts the algorithm.
deba@948
  1200
    ///
deba@948
  1201
    /// \pre \ref init() must be called before using this function.
deba@948
  1202
    void start() {
deba@948
  1203
      enum OpType {
deba@948
  1204
        D1, D2, D3
deba@948
  1205
      };
deba@948
  1206
deba@948
  1207
      int unmatched = _node_num;
deba@948
  1208
      while (unmatched > 0) {
deba@948
  1209
        Value d1 = !_delta1->empty() ?
deba@948
  1210
          _delta1->prio() : std::numeric_limits<Value>::max();
deba@948
  1211
deba@948
  1212
        Value d2 = !_delta2->empty() ?
deba@948
  1213
          _delta2->prio() : std::numeric_limits<Value>::max();
deba@948
  1214
deba@948
  1215
        Value d3 = !_delta3->empty() ?
deba@948
  1216
          _delta3->prio() : std::numeric_limits<Value>::max();
deba@948
  1217
deba@948
  1218
        _delta_sum = d3; OpType ot = D3;
deba@948
  1219
        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
deba@948
  1220
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
deba@948
  1221
deba@948
  1222
        switch (ot) {
deba@948
  1223
        case D1:
deba@948
  1224
          {
deba@948
  1225
            Node n = _delta1->top();
deba@948
  1226
            unmatchNode(n);
deba@948
  1227
            --unmatched;
deba@948
  1228
          }
deba@948
  1229
          break;
deba@948
  1230
        case D2:
deba@948
  1231
          {
deba@948
  1232
            Node n = _delta2->top();
deba@948
  1233
            Arc a = (*_pred)[n];
deba@948
  1234
            if ((*_matching)[n] == INVALID) {
deba@948
  1235
              augmentOnArc(a);
deba@948
  1236
              --unmatched;
deba@948
  1237
            } else {
deba@948
  1238
              Node v = _graph.target((*_matching)[n]);
deba@948
  1239
              if ((*_matching)[n] !=
deba@948
  1240
                  _graph.oppositeArc((*_matching)[v])) {
deba@948
  1241
                extractCycle(a);
deba@948
  1242
                --unmatched;
deba@948
  1243
              } else {
deba@948
  1244
                extendOnArc(a);
deba@948
  1245
              }
deba@948
  1246
            }
deba@948
  1247
          } break;
deba@948
  1248
        case D3:
deba@948
  1249
          {
deba@948
  1250
            Edge e = _delta3->top();
deba@948
  1251
deba@948
  1252
            Node left = _graph.u(e);
deba@948
  1253
            Node right = _graph.v(e);
deba@948
  1254
deba@948
  1255
            int left_tree = _tree_set->find(left);
deba@948
  1256
            int right_tree = _tree_set->find(right);
deba@948
  1257
deba@948
  1258
            if (left_tree == right_tree) {
deba@948
  1259
              cycleOnEdge(e, left_tree);
deba@948
  1260
              --unmatched;
deba@948
  1261
            } else {
deba@948
  1262
              augmentOnEdge(e);
deba@948
  1263
              unmatched -= 2;
deba@948
  1264
            }
deba@948
  1265
          } break;
deba@948
  1266
        }
deba@948
  1267
      }
deba@948
  1268
    }
deba@948
  1269
deba@948
  1270
    /// \brief Run the algorithm.
deba@948
  1271
    ///
deba@950
  1272
    /// This method runs the \c %MaxWeightedFractionalMatching algorithm.
deba@948
  1273
    ///
deba@948
  1274
    /// \note mwfm.run() is just a shortcut of the following code.
deba@948
  1275
    /// \code
deba@948
  1276
    ///   mwfm.init();
deba@948
  1277
    ///   mwfm.start();
deba@948
  1278
    /// \endcode
deba@948
  1279
    void run() {
deba@948
  1280
      init();
deba@948
  1281
      start();
deba@948
  1282
    }
deba@948
  1283
deba@948
  1284
    /// @}
deba@948
  1285
deba@948
  1286
    /// \name Primal Solution
deba@948
  1287
    /// Functions to get the primal solution, i.e. the maximum weighted
deba@948
  1288
    /// matching.\n
deba@948
  1289
    /// Either \ref run() or \ref start() function should be called before
deba@948
  1290
    /// using them.
deba@948
  1291
deba@948
  1292
    /// @{
deba@948
  1293
deba@948
  1294
    /// \brief Return the weight of the matching.
deba@948
  1295
    ///
deba@948
  1296
    /// This function returns the weight of the found matching. This
deba@948
  1297
    /// value is scaled by \ref primalScale "primal scale".
deba@948
  1298
    ///
deba@948
  1299
    /// \pre Either run() or start() must be called before using this function.
deba@948
  1300
    Value matchingWeight() const {
deba@948
  1301
      Value sum = 0;
deba@948
  1302
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
  1303
        if ((*_matching)[n] != INVALID) {
deba@948
  1304
          sum += _weight[(*_matching)[n]];
deba@948
  1305
        }
deba@948
  1306
      }
deba@948
  1307
      return sum * primalScale / 2;
deba@948
  1308
    }
deba@948
  1309
deba@948
  1310
    /// \brief Return the number of covered nodes in the matching.
deba@948
  1311
    ///
deba@948
  1312
    /// This function returns the number of covered nodes in the matching.
deba@948
  1313
    ///
deba@948
  1314
    /// \pre Either run() or start() must be called before using this function.
deba@948
  1315
    int matchingSize() const {
deba@948
  1316
      int num = 0;
deba@948
  1317
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
  1318
        if ((*_matching)[n] != INVALID) {
deba@948
  1319
          ++num;
deba@948
  1320
        }
deba@948
  1321
      }
deba@948
  1322
      return num;
deba@948
  1323
    }
deba@948
  1324
deba@948
  1325
    /// \brief Return \c true if the given edge is in the matching.
deba@948
  1326
    ///
deba@948
  1327
    /// This function returns \c true if the given edge is in the
deba@948
  1328
    /// found matching. The result is scaled by \ref primalScale
deba@948
  1329
    /// "primal scale".
deba@948
  1330
    ///
deba@948
  1331
    /// \pre Either run() or start() must be called before using this function.
deba@948
  1332
    Value matching(const Edge& edge) const {
deba@948
  1333
      return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
deba@948
  1334
        * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
deba@948
  1335
        * primalScale / 2;
deba@948
  1336
    }
deba@948
  1337
deba@948
  1338
    /// \brief Return the fractional matching arc (or edge) incident
deba@948
  1339
    /// to the given node.
deba@948
  1340
    ///
deba@948
  1341
    /// This function returns one of the fractional matching arc (or
deba@948
  1342
    /// edge) incident to the given node in the found matching or \c
deba@948
  1343
    /// INVALID if the node is not covered by the matching or if the
deba@948
  1344
    /// node is on an odd length cycle then it is the successor edge
deba@948
  1345
    /// on the cycle.
deba@948
  1346
    ///
deba@948
  1347
    /// \pre Either run() or start() must be called before using this function.
deba@948
  1348
    Arc matching(const Node& node) const {
deba@948
  1349
      return (*_matching)[node];
deba@948
  1350
    }
deba@948
  1351
deba@948
  1352
    /// \brief Return a const reference to the matching map.
deba@948
  1353
    ///
deba@948
  1354
    /// This function returns a const reference to a node map that stores
deba@948
  1355
    /// the matching arc (or edge) incident to each node.
deba@948
  1356
    const MatchingMap& matchingMap() const {
deba@948
  1357
      return *_matching;
deba@948
  1358
    }
deba@948
  1359
deba@948
  1360
    /// @}
deba@948
  1361
deba@948
  1362
    /// \name Dual Solution
deba@948
  1363
    /// Functions to get the dual solution.\n
deba@948
  1364
    /// Either \ref run() or \ref start() function should be called before
deba@948
  1365
    /// using them.
deba@948
  1366
deba@948
  1367
    /// @{
deba@948
  1368
deba@948
  1369
    /// \brief Return the value of the dual solution.
deba@948
  1370
    ///
deba@948
  1371
    /// This function returns the value of the dual solution.
deba@948
  1372
    /// It should be equal to the primal value scaled by \ref dualScale
deba@948
  1373
    /// "dual scale".
deba@948
  1374
    ///
deba@948
  1375
    /// \pre Either run() or start() must be called before using this function.
deba@948
  1376
    Value dualValue() const {
deba@948
  1377
      Value sum = 0;
deba@948
  1378
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
  1379
        sum += nodeValue(n);
deba@948
  1380
      }
deba@948
  1381
      return sum;
deba@948
  1382
    }
deba@948
  1383
deba@948
  1384
    /// \brief Return the dual value (potential) of the given node.
deba@948
  1385
    ///
deba@948
  1386
    /// This function returns the dual value (potential) of the given node.
deba@948
  1387
    ///
deba@948
  1388
    /// \pre Either run() or start() must be called before using this function.
deba@948
  1389
    Value nodeValue(const Node& n) const {
deba@948
  1390
      return (*_node_potential)[n];
deba@948
  1391
    }
deba@948
  1392
deba@948
  1393
    /// @}
deba@948
  1394
deba@948
  1395
  };
deba@948
  1396
deba@948
  1397
  /// \ingroup matching
deba@948
  1398
  ///
deba@948
  1399
  /// \brief Weighted fractional perfect matching in general graphs
deba@948
  1400
  ///
deba@948
  1401
  /// This class provides an efficient implementation of fractional
deba@950
  1402
  /// matching algorithm. The implementation uses priority queues and
deba@950
  1403
  /// provides \f$O(nm\log n)\f$ time complexity.
deba@948
  1404
  ///
deba@948
  1405
  /// The maximum weighted fractional perfect matching is a relaxation
deba@948
  1406
  /// of the maximum weighted perfect matching problem where the odd
deba@948
  1407
  /// set constraints are omitted.
deba@948
  1408
  /// It can be formulated with the following linear program.
deba@948
  1409
  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
deba@948
  1410
  /// \f[x_e \ge 0\quad \forall e\in E\f]
deba@948
  1411
  /// \f[\max \sum_{e\in E}x_ew_e\f]
deba@948
  1412
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
deba@948
  1413
  /// \f$X\f$. The result must be the union of a matching with one
deba@948
  1414
  /// value edges and a set of odd length cycles with half value edges.
deba@948
  1415
  ///
deba@948
  1416
  /// The algorithm calculates an optimal fractional matching and a
deba@948
  1417
  /// proof of the optimality. The solution of the dual problem can be
deba@948
  1418
  /// used to check the result of the algorithm. The dual linear
deba@948
  1419
  /// problem is the following.
deba@948
  1420
  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
deba@950
  1421
  /// \f[\min \sum_{u \in V}y_u \f]
deba@948
  1422
  ///
deba@948
  1423
  /// The algorithm can be executed with the run() function.
deba@948
  1424
  /// After it the matching (the primal solution) and the dual solution
deba@948
  1425
  /// can be obtained using the query functions.
deba@948
  1426
deba@948
  1427
  /// If the value type is integer, then the primal and the dual
deba@948
  1428
  /// solutions are multiplied by
deba@950
  1429
  /// \ref MaxWeightedPerfectFractionalMatching::primalScale "2" and
deba@950
  1430
  /// \ref MaxWeightedPerfectFractionalMatching::dualScale "4" respectively.
deba@948
  1431
  ///
deba@948
  1432
  /// \tparam GR The undirected graph type the algorithm runs on.
deba@948
  1433
  /// \tparam WM The type edge weight map. The default type is
deba@948
  1434
  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
deba@948
  1435
#ifdef DOXYGEN
deba@948
  1436
  template <typename GR, typename WM>
deba@948
  1437
#else
deba@948
  1438
  template <typename GR,
deba@948
  1439
            typename WM = typename GR::template EdgeMap<int> >
deba@948
  1440
#endif
deba@948
  1441
  class MaxWeightedPerfectFractionalMatching {
deba@948
  1442
  public:
deba@948
  1443
deba@948
  1444
    /// The graph type of the algorithm
deba@948
  1445
    typedef GR Graph;
deba@948
  1446
    /// The type of the edge weight map
deba@948
  1447
    typedef WM WeightMap;
deba@948
  1448
    /// The value type of the edge weights
deba@948
  1449
    typedef typename WeightMap::Value Value;
deba@948
  1450
deba@948
  1451
    /// The type of the matching map
deba@948
  1452
    typedef typename Graph::template NodeMap<typename Graph::Arc>
deba@948
  1453
    MatchingMap;
deba@948
  1454
deba@948
  1455
    /// \brief Scaling factor for primal solution
deba@948
  1456
    ///
deba@948
  1457
    /// Scaling factor for primal solution. It is equal to 2 or 1
deba@948
  1458
    /// according to the value type.
deba@948
  1459
    static const int primalScale =
deba@948
  1460
      std::numeric_limits<Value>::is_integer ? 2 : 1;
deba@948
  1461
deba@948
  1462
    /// \brief Scaling factor for dual solution
deba@948
  1463
    ///
deba@948
  1464
    /// Scaling factor for dual solution. It is equal to 4 or 1
deba@948
  1465
    /// according to the value type.
deba@948
  1466
    static const int dualScale =
deba@948
  1467
      std::numeric_limits<Value>::is_integer ? 4 : 1;
deba@948
  1468
deba@948
  1469
  private:
deba@948
  1470
deba@948
  1471
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
deba@948
  1472
deba@948
  1473
    typedef typename Graph::template NodeMap<Value> NodePotential;
deba@948
  1474
deba@948
  1475
    const Graph& _graph;
deba@948
  1476
    const WeightMap& _weight;
deba@948
  1477
deba@948
  1478
    MatchingMap* _matching;
deba@948
  1479
    NodePotential* _node_potential;
deba@948
  1480
deba@948
  1481
    int _node_num;
deba@948
  1482
    bool _allow_loops;
deba@948
  1483
deba@948
  1484
    enum Status {
deba@948
  1485
      EVEN = -1, MATCHED = 0, ODD = 1
deba@948
  1486
    };
deba@948
  1487
deba@948
  1488
    typedef typename Graph::template NodeMap<Status> StatusMap;
deba@948
  1489
    StatusMap* _status;
deba@948
  1490
deba@948
  1491
    typedef typename Graph::template NodeMap<Arc> PredMap;
deba@948
  1492
    PredMap* _pred;
deba@948
  1493
deba@948
  1494
    typedef ExtendFindEnum<IntNodeMap> TreeSet;
deba@948
  1495
deba@948
  1496
    IntNodeMap *_tree_set_index;
deba@948
  1497
    TreeSet *_tree_set;
deba@948
  1498
deba@948
  1499
    IntNodeMap *_delta2_index;
deba@948
  1500
    BinHeap<Value, IntNodeMap> *_delta2;
deba@948
  1501
deba@948
  1502
    IntEdgeMap *_delta3_index;
deba@948
  1503
    BinHeap<Value, IntEdgeMap> *_delta3;
deba@948
  1504
deba@948
  1505
    Value _delta_sum;
deba@948
  1506
deba@948
  1507
    void createStructures() {
deba@948
  1508
      _node_num = countNodes(_graph);
deba@948
  1509
deba@948
  1510
      if (!_matching) {
deba@948
  1511
        _matching = new MatchingMap(_graph);
deba@948
  1512
      }
deba@948
  1513
      if (!_node_potential) {
deba@948
  1514
        _node_potential = new NodePotential(_graph);
deba@948
  1515
      }
deba@948
  1516
      if (!_status) {
deba@948
  1517
        _status = new StatusMap(_graph);
deba@948
  1518
      }
deba@948
  1519
      if (!_pred) {
deba@948
  1520
        _pred = new PredMap(_graph);
deba@948
  1521
      }
deba@948
  1522
      if (!_tree_set) {
deba@948
  1523
        _tree_set_index = new IntNodeMap(_graph);
deba@948
  1524
        _tree_set = new TreeSet(*_tree_set_index);
deba@948
  1525
      }
deba@948
  1526
      if (!_delta2) {
deba@948
  1527
        _delta2_index = new IntNodeMap(_graph);
deba@948
  1528
        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
deba@948
  1529
      }
deba@948
  1530
      if (!_delta3) {
deba@948
  1531
        _delta3_index = new IntEdgeMap(_graph);
deba@948
  1532
        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
deba@948
  1533
      }
deba@948
  1534
    }
deba@948
  1535
deba@948
  1536
    void destroyStructures() {
deba@948
  1537
      if (_matching) {
deba@948
  1538
        delete _matching;
deba@948
  1539
      }
deba@948
  1540
      if (_node_potential) {
deba@948
  1541
        delete _node_potential;
deba@948
  1542
      }
deba@948
  1543
      if (_status) {
deba@948
  1544
        delete _status;
deba@948
  1545
      }
deba@948
  1546
      if (_pred) {
deba@948
  1547
        delete _pred;
deba@948
  1548
      }
deba@948
  1549
      if (_tree_set) {
deba@948
  1550
        delete _tree_set_index;
deba@948
  1551
        delete _tree_set;
deba@948
  1552
      }
deba@948
  1553
      if (_delta2) {
deba@948
  1554
        delete _delta2_index;
deba@948
  1555
        delete _delta2;
deba@948
  1556
      }
deba@948
  1557
      if (_delta3) {
deba@948
  1558
        delete _delta3_index;
deba@948
  1559
        delete _delta3;
deba@948
  1560
      }
deba@948
  1561
    }
deba@948
  1562
deba@948
  1563
    void matchedToEven(Node node, int tree) {
deba@948
  1564
      _tree_set->insert(node, tree);
deba@948
  1565
      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
deba@948
  1566
deba@948
  1567
      if (_delta2->state(node) == _delta2->IN_HEAP) {
deba@948
  1568
        _delta2->erase(node);
deba@948
  1569
      }
deba@948
  1570
deba@948
  1571
      for (InArcIt a(_graph, node); a != INVALID; ++a) {
deba@948
  1572
        Node v = _graph.source(a);
deba@948
  1573
        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
deba@948
  1574
          dualScale * _weight[a];
deba@948
  1575
        if (node == v) {
deba@948
  1576
          if (_allow_loops && _graph.direction(a)) {
deba@948
  1577
            _delta3->push(a, rw / 2);
deba@948
  1578
          }
deba@948
  1579
        } else if ((*_status)[v] == EVEN) {
deba@948
  1580
          _delta3->push(a, rw / 2);
deba@948
  1581
        } else if ((*_status)[v] == MATCHED) {
deba@948
  1582
          if (_delta2->state(v) != _delta2->IN_HEAP) {
deba@948
  1583
            _pred->set(v, a);
deba@948
  1584
            _delta2->push(v, rw);
deba@948
  1585
          } else if ((*_delta2)[v] > rw) {
deba@948
  1586
            _pred->set(v, a);
deba@948
  1587
            _delta2->decrease(v, rw);
deba@948
  1588
          }
deba@948
  1589
        }
deba@948
  1590
      }
deba@948
  1591
    }
deba@948
  1592
deba@948
  1593
    void matchedToOdd(Node node, int tree) {
deba@948
  1594
      _tree_set->insert(node, tree);
deba@948
  1595
      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
deba@948
  1596
deba@948
  1597
      if (_delta2->state(node) == _delta2->IN_HEAP) {
deba@948
  1598
        _delta2->erase(node);
deba@948
  1599
      }
deba@948
  1600
    }
deba@948
  1601
deba@948
  1602
    void evenToMatched(Node node, int tree) {
deba@948
  1603
      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
deba@948
  1604
      Arc min = INVALID;
deba@948
  1605
      Value minrw = std::numeric_limits<Value>::max();
deba@948
  1606
      for (InArcIt a(_graph, node); a != INVALID; ++a) {
deba@948
  1607
        Node v = _graph.source(a);
deba@948
  1608
        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
deba@948
  1609
          dualScale * _weight[a];
deba@948
  1610
deba@948
  1611
        if (node == v) {
deba@948
  1612
          if (_allow_loops && _graph.direction(a)) {
deba@948
  1613
            _delta3->erase(a);
deba@948
  1614
          }
deba@948
  1615
        } else if ((*_status)[v] == EVEN) {
deba@948
  1616
          _delta3->erase(a);
deba@948
  1617
          if (minrw > rw) {
deba@948
  1618
            min = _graph.oppositeArc(a);
deba@948
  1619
            minrw = rw;
deba@948
  1620
          }
deba@948
  1621
        } else if ((*_status)[v]  == MATCHED) {
deba@948
  1622
          if ((*_pred)[v] == a) {
deba@948
  1623
            Arc mina = INVALID;
deba@948
  1624
            Value minrwa = std::numeric_limits<Value>::max();
deba@948
  1625
            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
deba@948
  1626
              Node va = _graph.target(aa);
deba@948
  1627
              if ((*_status)[va] != EVEN ||
deba@948
  1628
                  _tree_set->find(va) == tree) continue;
deba@948
  1629
              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
deba@948
  1630
                dualScale * _weight[aa];
deba@948
  1631
              if (minrwa > rwa) {
deba@948
  1632
                minrwa = rwa;
deba@948
  1633
                mina = aa;
deba@948
  1634
              }
deba@948
  1635
            }
deba@948
  1636
            if (mina != INVALID) {
deba@948
  1637
              _pred->set(v, mina);
deba@948
  1638
              _delta2->increase(v, minrwa);
deba@948
  1639
            } else {
deba@948
  1640
              _pred->set(v, INVALID);
deba@948
  1641
              _delta2->erase(v);
deba@948
  1642
            }
deba@948
  1643
          }
deba@948
  1644
        }
deba@948
  1645
      }
deba@948
  1646
      if (min != INVALID) {
deba@948
  1647
        _pred->set(node, min);
deba@948
  1648
        _delta2->push(node, minrw);
deba@948
  1649
      } else {
deba@948
  1650
        _pred->set(node, INVALID);
deba@948
  1651
      }
deba@948
  1652
    }
deba@948
  1653
deba@948
  1654
    void oddToMatched(Node node) {
deba@948
  1655
      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
deba@948
  1656
      Arc min = INVALID;
deba@948
  1657
      Value minrw = std::numeric_limits<Value>::max();
deba@948
  1658
      for (InArcIt a(_graph, node); a != INVALID; ++a) {
deba@948
  1659
        Node v = _graph.source(a);
deba@948
  1660
        if ((*_status)[v] != EVEN) continue;
deba@948
  1661
        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
deba@948
  1662
          dualScale * _weight[a];
deba@948
  1663
deba@948
  1664
        if (minrw > rw) {
deba@948
  1665
          min = _graph.oppositeArc(a);
deba@948
  1666
          minrw = rw;
deba@948
  1667
        }
deba@948
  1668
      }
deba@948
  1669
      if (min != INVALID) {
deba@948
  1670
        _pred->set(node, min);
deba@948
  1671
        _delta2->push(node, minrw);
deba@948
  1672
      } else {
deba@948
  1673
        _pred->set(node, INVALID);
deba@948
  1674
      }
deba@948
  1675
    }
deba@948
  1676
deba@948
  1677
    void alternatePath(Node even, int tree) {
deba@948
  1678
      Node odd;
deba@948
  1679
deba@948
  1680
      _status->set(even, MATCHED);
deba@948
  1681
      evenToMatched(even, tree);
deba@948
  1682
deba@948
  1683
      Arc prev = (*_matching)[even];
deba@948
  1684
      while (prev != INVALID) {
deba@948
  1685
        odd = _graph.target(prev);
deba@948
  1686
        even = _graph.target((*_pred)[odd]);
deba@948
  1687
        _matching->set(odd, (*_pred)[odd]);
deba@948
  1688
        _status->set(odd, MATCHED);
deba@948
  1689
        oddToMatched(odd);
deba@948
  1690
deba@948
  1691
        prev = (*_matching)[even];
deba@948
  1692
        _status->set(even, MATCHED);
deba@948
  1693
        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
deba@948
  1694
        evenToMatched(even, tree);
deba@948
  1695
      }
deba@948
  1696
    }
deba@948
  1697
deba@948
  1698
    void destroyTree(int tree) {
deba@948
  1699
      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
deba@948
  1700
        if ((*_status)[n] == EVEN) {
deba@948
  1701
          _status->set(n, MATCHED);
deba@948
  1702
          evenToMatched(n, tree);
deba@948
  1703
        } else if ((*_status)[n] == ODD) {
deba@948
  1704
          _status->set(n, MATCHED);
deba@948
  1705
          oddToMatched(n);
deba@948
  1706
        }
deba@948
  1707
      }
deba@948
  1708
      _tree_set->eraseClass(tree);
deba@948
  1709
    }
deba@948
  1710
deba@948
  1711
    void augmentOnEdge(const Edge& edge) {
deba@948
  1712
      Node left = _graph.u(edge);
deba@948
  1713
      int left_tree = _tree_set->find(left);
deba@948
  1714
deba@948
  1715
      alternatePath(left, left_tree);
deba@948
  1716
      destroyTree(left_tree);
deba@948
  1717
      _matching->set(left, _graph.direct(edge, true));
deba@948
  1718
deba@948
  1719
      Node right = _graph.v(edge);
deba@948
  1720
      int right_tree = _tree_set->find(right);
deba@948
  1721
deba@948
  1722
      alternatePath(right, right_tree);
deba@948
  1723
      destroyTree(right_tree);
deba@948
  1724
      _matching->set(right, _graph.direct(edge, false));
deba@948
  1725
    }
deba@948
  1726
deba@948
  1727
    void augmentOnArc(const Arc& arc) {
deba@948
  1728
      Node left = _graph.source(arc);
deba@948
  1729
      _status->set(left, MATCHED);
deba@948
  1730
      _matching->set(left, arc);
deba@948
  1731
      _pred->set(left, arc);
deba@948
  1732
deba@948
  1733
      Node right = _graph.target(arc);
deba@948
  1734
      int right_tree = _tree_set->find(right);
deba@948
  1735
deba@948
  1736
      alternatePath(right, right_tree);
deba@948
  1737
      destroyTree(right_tree);
deba@948
  1738
      _matching->set(right, _graph.oppositeArc(arc));
deba@948
  1739
    }
deba@948
  1740
deba@948
  1741
    void extendOnArc(const Arc& arc) {
deba@948
  1742
      Node base = _graph.target(arc);
deba@948
  1743
      int tree = _tree_set->find(base);
deba@948
  1744
deba@948
  1745
      Node odd = _graph.source(arc);
deba@948
  1746
      _tree_set->insert(odd, tree);
deba@948
  1747
      _status->set(odd, ODD);
deba@948
  1748
      matchedToOdd(odd, tree);
deba@948
  1749
      _pred->set(odd, arc);
deba@948
  1750
deba@948
  1751
      Node even = _graph.target((*_matching)[odd]);
deba@948
  1752
      _tree_set->insert(even, tree);
deba@948
  1753
      _status->set(even, EVEN);
deba@948
  1754
      matchedToEven(even, tree);
deba@948
  1755
    }
deba@948
  1756
deba@948
  1757
    void cycleOnEdge(const Edge& edge, int tree) {
deba@948
  1758
      Node nca = INVALID;
deba@948
  1759
      std::vector<Node> left_path, right_path;
deba@948
  1760
deba@948
  1761
      {
deba@948
  1762
        std::set<Node> left_set, right_set;
deba@948
  1763
        Node left = _graph.u(edge);
deba@948
  1764
        left_path.push_back(left);
deba@948
  1765
        left_set.insert(left);
deba@948
  1766
deba@948
  1767
        Node right = _graph.v(edge);
deba@948
  1768
        right_path.push_back(right);
deba@948
  1769
        right_set.insert(right);
deba@948
  1770
deba@948
  1771
        while (true) {
deba@948
  1772
deba@948
  1773
          if (left_set.find(right) != left_set.end()) {
deba@948
  1774
            nca = right;
deba@948
  1775
            break;
deba@948
  1776
          }
deba@948
  1777
deba@948
  1778
          if ((*_matching)[left] == INVALID) break;
deba@948
  1779
deba@948
  1780
          left = _graph.target((*_matching)[left]);
deba@948
  1781
          left_path.push_back(left);
deba@948
  1782
          left = _graph.target((*_pred)[left]);
deba@948
  1783
          left_path.push_back(left);
deba@948
  1784
deba@948
  1785
          left_set.insert(left);
deba@948
  1786
deba@948
  1787
          if (right_set.find(left) != right_set.end()) {
deba@948
  1788
            nca = left;
deba@948
  1789
            break;
deba@948
  1790
          }
deba@948
  1791
deba@948
  1792
          if ((*_matching)[right] == INVALID) break;
deba@948
  1793
deba@948
  1794
          right = _graph.target((*_matching)[right]);
deba@948
  1795
          right_path.push_back(right);
deba@948
  1796
          right = _graph.target((*_pred)[right]);
deba@948
  1797
          right_path.push_back(right);
deba@948
  1798
deba@948
  1799
          right_set.insert(right);
deba@948
  1800
deba@948
  1801
        }
deba@948
  1802
deba@948
  1803
        if (nca == INVALID) {
deba@948
  1804
          if ((*_matching)[left] == INVALID) {
deba@948
  1805
            nca = right;
deba@948
  1806
            while (left_set.find(nca) == left_set.end()) {
deba@948
  1807
              nca = _graph.target((*_matching)[nca]);
deba@948
  1808
              right_path.push_back(nca);
deba@948
  1809
              nca = _graph.target((*_pred)[nca]);
deba@948
  1810
              right_path.push_back(nca);
deba@948
  1811
            }
deba@948
  1812
          } else {
deba@948
  1813
            nca = left;
deba@948
  1814
            while (right_set.find(nca) == right_set.end()) {
deba@948
  1815
              nca = _graph.target((*_matching)[nca]);
deba@948
  1816
              left_path.push_back(nca);
deba@948
  1817
              nca = _graph.target((*_pred)[nca]);
deba@948
  1818
              left_path.push_back(nca);
deba@948
  1819
            }
deba@948
  1820
          }
deba@948
  1821
        }
deba@948
  1822
      }
deba@948
  1823
deba@948
  1824
      alternatePath(nca, tree);
deba@948
  1825
      Arc prev;
deba@948
  1826
deba@948
  1827
      prev = _graph.direct(edge, true);
deba@948
  1828
      for (int i = 0; left_path[i] != nca; i += 2) {
deba@948
  1829
        _matching->set(left_path[i], prev);
deba@948
  1830
        _status->set(left_path[i], MATCHED);
deba@948
  1831
        evenToMatched(left_path[i], tree);
deba@948
  1832
deba@948
  1833
        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
deba@948
  1834
        _status->set(left_path[i + 1], MATCHED);
deba@948
  1835
        oddToMatched(left_path[i + 1]);
deba@948
  1836
      }
deba@948
  1837
      _matching->set(nca, prev);
deba@948
  1838
deba@948
  1839
      for (int i = 0; right_path[i] != nca; i += 2) {
deba@948
  1840
        _status->set(right_path[i], MATCHED);
deba@948
  1841
        evenToMatched(right_path[i], tree);
deba@948
  1842
deba@948
  1843
        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
deba@948
  1844
        _status->set(right_path[i + 1], MATCHED);
deba@948
  1845
        oddToMatched(right_path[i + 1]);
deba@948
  1846
      }
deba@948
  1847
deba@948
  1848
      destroyTree(tree);
deba@948
  1849
    }
deba@948
  1850
deba@948
  1851
    void extractCycle(const Arc &arc) {
deba@948
  1852
      Node left = _graph.source(arc);
deba@948
  1853
      Node odd = _graph.target((*_matching)[left]);
deba@948
  1854
      Arc prev;
deba@948
  1855
      while (odd != left) {
deba@948
  1856
        Node even = _graph.target((*_matching)[odd]);
deba@948
  1857
        prev = (*_matching)[odd];
deba@948
  1858
        odd = _graph.target((*_matching)[even]);
deba@948
  1859
        _matching->set(even, _graph.oppositeArc(prev));
deba@948
  1860
      }
deba@948
  1861
      _matching->set(left, arc);
deba@948
  1862
deba@948
  1863
      Node right = _graph.target(arc);
deba@948
  1864
      int right_tree = _tree_set->find(right);
deba@948
  1865
      alternatePath(right, right_tree);
deba@948
  1866
      destroyTree(right_tree);
deba@948
  1867
      _matching->set(right, _graph.oppositeArc(arc));
deba@948
  1868
    }
deba@948
  1869
deba@948
  1870
  public:
deba@948
  1871
deba@948
  1872
    /// \brief Constructor
deba@948
  1873
    ///
deba@948
  1874
    /// Constructor.
deba@948
  1875
    MaxWeightedPerfectFractionalMatching(const Graph& graph,
deba@948
  1876
                                         const WeightMap& weight,
deba@948
  1877
                                         bool allow_loops = true)
deba@948
  1878
      : _graph(graph), _weight(weight), _matching(0),
deba@948
  1879
      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
deba@948
  1880
      _status(0),  _pred(0),
deba@948
  1881
      _tree_set_index(0), _tree_set(0),
deba@948
  1882
deba@948
  1883
      _delta2_index(0), _delta2(0),
deba@948
  1884
      _delta3_index(0), _delta3(0),
deba@948
  1885
deba@948
  1886
      _delta_sum() {}
deba@948
  1887
deba@948
  1888
    ~MaxWeightedPerfectFractionalMatching() {
deba@948
  1889
      destroyStructures();
deba@948
  1890
    }
deba@948
  1891
deba@948
  1892
    /// \name Execution Control
deba@948
  1893
    /// The simplest way to execute the algorithm is to use the
deba@948
  1894
    /// \ref run() member function.
deba@948
  1895
deba@948
  1896
    ///@{
deba@948
  1897
deba@948
  1898
    /// \brief Initialize the algorithm
deba@948
  1899
    ///
deba@948
  1900
    /// This function initializes the algorithm.
deba@948
  1901
    void init() {
deba@948
  1902
      createStructures();
deba@948
  1903
deba@948
  1904
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
  1905
        (*_delta2_index)[n] = _delta2->PRE_HEAP;
deba@948
  1906
      }
deba@948
  1907
      for (EdgeIt e(_graph); e != INVALID; ++e) {
deba@948
  1908
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
deba@948
  1909
      }
deba@948
  1910
deba@948
  1911
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
  1912
        Value max = - std::numeric_limits<Value>::max();
deba@948
  1913
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
deba@948
  1914
          if (_graph.target(e) == n && !_allow_loops) continue;
deba@948
  1915
          if ((dualScale * _weight[e]) / 2 > max) {
deba@948
  1916
            max = (dualScale * _weight[e]) / 2;
deba@948
  1917
          }
deba@948
  1918
        }
deba@948
  1919
        _node_potential->set(n, max);
deba@948
  1920
deba@948
  1921
        _tree_set->insert(n);
deba@948
  1922
deba@948
  1923
        _matching->set(n, INVALID);
deba@948
  1924
        _status->set(n, EVEN);
deba@948
  1925
      }
deba@948
  1926
deba@948
  1927
      for (EdgeIt e(_graph); e != INVALID; ++e) {
deba@948
  1928
        Node left = _graph.u(e);
deba@948
  1929
        Node right = _graph.v(e);
deba@948
  1930
        if (left == right && !_allow_loops) continue;
deba@948
  1931
        _delta3->push(e, ((*_node_potential)[left] +
deba@948
  1932
                          (*_node_potential)[right] -
deba@948
  1933
                          dualScale * _weight[e]) / 2);
deba@948
  1934
      }
deba@948
  1935
    }
deba@948
  1936
deba@948
  1937
    /// \brief Start the algorithm
deba@948
  1938
    ///
deba@948
  1939
    /// This function starts the algorithm.
deba@948
  1940
    ///
deba@948
  1941
    /// \pre \ref init() must be called before using this function.
deba@948
  1942
    bool start() {
deba@948
  1943
      enum OpType {
deba@948
  1944
        D2, D3
deba@948
  1945
      };
deba@948
  1946
deba@948
  1947
      int unmatched = _node_num;
deba@948
  1948
      while (unmatched > 0) {
deba@948
  1949
        Value d2 = !_delta2->empty() ?
deba@948
  1950
          _delta2->prio() : std::numeric_limits<Value>::max();
deba@948
  1951
deba@948
  1952
        Value d3 = !_delta3->empty() ?
deba@948
  1953
          _delta3->prio() : std::numeric_limits<Value>::max();
deba@948
  1954
deba@948
  1955
        _delta_sum = d3; OpType ot = D3;
deba@948
  1956
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
deba@948
  1957
deba@948
  1958
        if (_delta_sum == std::numeric_limits<Value>::max()) {
deba@948
  1959
          return false;
deba@948
  1960
        }
deba@948
  1961
deba@948
  1962
        switch (ot) {
deba@948
  1963
        case D2:
deba@948
  1964
          {
deba@948
  1965
            Node n = _delta2->top();
deba@948
  1966
            Arc a = (*_pred)[n];
deba@948
  1967
            if ((*_matching)[n] == INVALID) {
deba@948
  1968
              augmentOnArc(a);
deba@948
  1969
              --unmatched;
deba@948
  1970
            } else {
deba@948
  1971
              Node v = _graph.target((*_matching)[n]);
deba@948
  1972
              if ((*_matching)[n] !=
deba@948
  1973
                  _graph.oppositeArc((*_matching)[v])) {
deba@948
  1974
                extractCycle(a);
deba@948
  1975
                --unmatched;
deba@948
  1976
              } else {
deba@948
  1977
                extendOnArc(a);
deba@948
  1978
              }
deba@948
  1979
            }
deba@948
  1980
          } break;
deba@948
  1981
        case D3:
deba@948
  1982
          {
deba@948
  1983
            Edge e = _delta3->top();
deba@948
  1984
deba@948
  1985
            Node left = _graph.u(e);
deba@948
  1986
            Node right = _graph.v(e);
deba@948
  1987
deba@948
  1988
            int left_tree = _tree_set->find(left);
deba@948
  1989
            int right_tree = _tree_set->find(right);
deba@948
  1990
deba@948
  1991
            if (left_tree == right_tree) {
deba@948
  1992
              cycleOnEdge(e, left_tree);
deba@948
  1993
              --unmatched;
deba@948
  1994
            } else {
deba@948
  1995
              augmentOnEdge(e);
deba@948
  1996
              unmatched -= 2;
deba@948
  1997
            }
deba@948
  1998
          } break;
deba@948
  1999
        }
deba@948
  2000
      }
deba@948
  2001
      return true;
deba@948
  2002
    }
deba@948
  2003
deba@948
  2004
    /// \brief Run the algorithm.
deba@948
  2005
    ///
deba@950
  2006
    /// This method runs the \c %MaxWeightedPerfectFractionalMatching 
deba@950
  2007
    /// algorithm.
deba@948
  2008
    ///
deba@948
  2009
    /// \note mwfm.run() is just a shortcut of the following code.
deba@948
  2010
    /// \code
deba@948
  2011
    ///   mwpfm.init();
deba@948
  2012
    ///   mwpfm.start();
deba@948
  2013
    /// \endcode
deba@948
  2014
    bool run() {
deba@948
  2015
      init();
deba@948
  2016
      return start();
deba@948
  2017
    }
deba@948
  2018
deba@948
  2019
    /// @}
deba@948
  2020
deba@948
  2021
    /// \name Primal Solution
deba@948
  2022
    /// Functions to get the primal solution, i.e. the maximum weighted
deba@948
  2023
    /// matching.\n
deba@948
  2024
    /// Either \ref run() or \ref start() function should be called before
deba@948
  2025
    /// using them.
deba@948
  2026
deba@948
  2027
    /// @{
deba@948
  2028
deba@948
  2029
    /// \brief Return the weight of the matching.
deba@948
  2030
    ///
deba@948
  2031
    /// This function returns the weight of the found matching. This
deba@948
  2032
    /// value is scaled by \ref primalScale "primal scale".
deba@948
  2033
    ///
deba@948
  2034
    /// \pre Either run() or start() must be called before using this function.
deba@948
  2035
    Value matchingWeight() const {
deba@948
  2036
      Value sum = 0;
deba@948
  2037
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
  2038
        if ((*_matching)[n] != INVALID) {
deba@948
  2039
          sum += _weight[(*_matching)[n]];
deba@948
  2040
        }
deba@948
  2041
      }
deba@948
  2042
      return sum * primalScale / 2;
deba@948
  2043
    }
deba@948
  2044
deba@948
  2045
    /// \brief Return the number of covered nodes in the matching.
deba@948
  2046
    ///
deba@948
  2047
    /// This function returns the number of covered nodes in the matching.
deba@948
  2048
    ///
deba@948
  2049
    /// \pre Either run() or start() must be called before using this function.
deba@948
  2050
    int matchingSize() const {
deba@948
  2051
      int num = 0;
deba@948
  2052
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
  2053
        if ((*_matching)[n] != INVALID) {
deba@948
  2054
          ++num;
deba@948
  2055
        }
deba@948
  2056
      }
deba@948
  2057
      return num;
deba@948
  2058
    }
deba@948
  2059
deba@948
  2060
    /// \brief Return \c true if the given edge is in the matching.
deba@948
  2061
    ///
deba@948
  2062
    /// This function returns \c true if the given edge is in the
deba@948
  2063
    /// found matching. The result is scaled by \ref primalScale
deba@948
  2064
    /// "primal scale".
deba@948
  2065
    ///
deba@948
  2066
    /// \pre Either run() or start() must be called before using this function.
deba@948
  2067
    Value matching(const Edge& edge) const {
deba@948
  2068
      return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
deba@948
  2069
        * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
deba@948
  2070
        * primalScale / 2;
deba@948
  2071
    }
deba@948
  2072
deba@948
  2073
    /// \brief Return the fractional matching arc (or edge) incident
deba@948
  2074
    /// to the given node.
deba@948
  2075
    ///
deba@948
  2076
    /// This function returns one of the fractional matching arc (or
deba@948
  2077
    /// edge) incident to the given node in the found matching or \c
deba@948
  2078
    /// INVALID if the node is not covered by the matching or if the
deba@948
  2079
    /// node is on an odd length cycle then it is the successor edge
deba@948
  2080
    /// on the cycle.
deba@948
  2081
    ///
deba@948
  2082
    /// \pre Either run() or start() must be called before using this function.
deba@948
  2083
    Arc matching(const Node& node) const {
deba@948
  2084
      return (*_matching)[node];
deba@948
  2085
    }
deba@948
  2086
deba@948
  2087
    /// \brief Return a const reference to the matching map.
deba@948
  2088
    ///
deba@948
  2089
    /// This function returns a const reference to a node map that stores
deba@948
  2090
    /// the matching arc (or edge) incident to each node.
deba@948
  2091
    const MatchingMap& matchingMap() const {
deba@948
  2092
      return *_matching;
deba@948
  2093
    }
deba@948
  2094
deba@948
  2095
    /// @}
deba@948
  2096
deba@948
  2097
    /// \name Dual Solution
deba@948
  2098
    /// Functions to get the dual solution.\n
deba@948
  2099
    /// Either \ref run() or \ref start() function should be called before
deba@948
  2100
    /// using them.
deba@948
  2101
deba@948
  2102
    /// @{
deba@948
  2103
deba@948
  2104
    /// \brief Return the value of the dual solution.
deba@948
  2105
    ///
deba@948
  2106
    /// This function returns the value of the dual solution.
deba@948
  2107
    /// It should be equal to the primal value scaled by \ref dualScale
deba@948
  2108
    /// "dual scale".
deba@948
  2109
    ///
deba@948
  2110
    /// \pre Either run() or start() must be called before using this function.
deba@948
  2111
    Value dualValue() const {
deba@948
  2112
      Value sum = 0;
deba@948
  2113
      for (NodeIt n(_graph); n != INVALID; ++n) {
deba@948
  2114
        sum += nodeValue(n);
deba@948
  2115
      }
deba@948
  2116
      return sum;
deba@948
  2117
    }
deba@948
  2118
deba@948
  2119
    /// \brief Return the dual value (potential) of the given node.
deba@948
  2120
    ///
deba@948
  2121
    /// This function returns the dual value (potential) of the given node.
deba@948
  2122
    ///
deba@948
  2123
    /// \pre Either run() or start() must be called before using this function.
deba@948
  2124
    Value nodeValue(const Node& n) const {
deba@948
  2125
      return (*_node_potential)[n];
deba@948
  2126
    }
deba@948
  2127
deba@948
  2128
    /// @}
deba@948
  2129
deba@948
  2130
  };
deba@948
  2131
deba@948
  2132
} //END OF NAMESPACE LEMON
deba@948
  2133
deba@948
  2134
#endif //LEMON_FRACTIONAL_MATCHING_H