COIN-OR::LEMON - Graph Library

Ticket #314: 636dadefe1e6.patch

File 636dadefe1e6.patch, 82.0 KB (added by Balazs Dezso, 11 years ago)

Fractional matching algorithms

  • doc/groups.dox

    # HG changeset patch
    # User Balazs Dezso <deba@inf.elte.hu>
    # Date 1253908296 -7200
    # Node ID 636dadefe1e6736f97988b05d6b4e7f95e4c5b9c
    # Parent  0513ccfea9670a94312cc2bde2ed8595f08361aa
    Add fractional matching algorithms (#314)
    
    diff -r 0513ccfea967 -r 636dadefe1e6 doc/groups.dox
    a b  
    349349also provide functions to query the minimum cut, which is the dual
    350350problem of maximum flow.
    351351
    352 \ref Circulation is a preflow push-relabel algorithm implemented directly 
     352\ref Circulation is a preflow push-relabel algorithm implemented directly
    353353for finding feasible circulations, which is a somewhat different problem,
    354354but it is strongly related to maximum flow.
    355355For more information, see \ref Circulation.
     
    470470- \ref MaxWeightedPerfectMatching
    471471  Edmond's blossom shrinking algorithm for calculating maximum weighted
    472472  perfect matching in general graphs.
     473- \ref MaxFractionalMatching Push-relabel algorithm for calculating
     474  maximum cardinality fractional matching in general graphs.
     475- \ref MaxWeightedFractionalMatching Augmenting path algorithm for calculating
     476  maximum weighted fractional matching in general graphs.
     477- \ref MaxWeightedPerfectFractionalMatching
     478  Augmenting path algorithm for calculating maximum weighted
     479  perfect fractional matching in general graphs.
    473480
    474481\image html bipartite_matching.png
    475482\image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth
  • lemon/Makefile.am

    diff -r 0513ccfea967 -r 636dadefe1e6 lemon/Makefile.am
    a b  
    8181        lemon/euler.h \
    8282        lemon/fib_heap.h \
    8383        lemon/fourary_heap.h \
     84        lemon/fractional_matching.h \
    8485        lemon/full_graph.h \
    8586        lemon/glpk.h \
    8687        lemon/gomory_hu.h \
  • new file lemon/fractional_matching.h

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

    diff -r 0513ccfea967 -r 636dadefe1e6 test/CMakeLists.txt
    a b  
    2121  edge_set_test
    2222  error_test
    2323  euler_test
     24  fractional_matching_test
    2425  gomory_hu_test
    2526  graph_copy_test
    2627  graph_test
  • test/Makefile.am

    diff -r 0513ccfea967 -r 636dadefe1e6 test/Makefile.am
    a b  
    1919        test/edge_set_test \
    2020        test/error_test \
    2121        test/euler_test \
     22        test/fractional_matching_test \
    2223        test/gomory_hu_test \
    2324        test/graph_copy_test \
    2425        test/graph_test \
     
    6566test_edge_set_test_SOURCES = test/edge_set_test.cc
    6667test_error_test_SOURCES = test/error_test.cc
    6768test_euler_test_SOURCES = test/euler_test.cc
     69test_fractional_matching_test_SOURCES = test/fractional_matching_test.cc
    6870test_gomory_hu_test_SOURCES = test/gomory_hu_test.cc
    6971test_graph_copy_test_SOURCES = test/graph_copy_test.cc
    7072test_graph_test_SOURCES = test/graph_test.cc
  • new file test/fractional_matching_test.cc

    diff -r 0513ccfea967 -r 636dadefe1e6 test/fractional_matching_test.cc
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2009
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <iostream>
     20#include <sstream>
     21#include <vector>
     22#include <queue>
     23#include <cstdlib>
     24
     25#include <lemon/fractional_matching.h>
     26#include <lemon/smart_graph.h>
     27#include <lemon/concepts/graph.h>
     28#include <lemon/concepts/maps.h>
     29#include <lemon/lgf_reader.h>
     30#include <lemon/math.h>
     31
     32#include "test_tools.h"
     33
     34using namespace std;
     35using namespace lemon;
     36
     37GRAPH_TYPEDEFS(SmartGraph);
     38
     39
     40const int lgfn = 4;
     41const std::string lgf[lgfn] = {
     42  "@nodes\n"
     43  "label\n"
     44  "0\n"
     45  "1\n"
     46  "2\n"
     47  "3\n"
     48  "4\n"
     49  "5\n"
     50  "6\n"
     51  "7\n"
     52  "@edges\n"
     53  "     label  weight\n"
     54  "7 4  0      984\n"
     55  "0 7  1      73\n"
     56  "7 1  2      204\n"
     57  "2 3  3      583\n"
     58  "2 7  4      565\n"
     59  "2 1  5      582\n"
     60  "0 4  6      551\n"
     61  "2 5  7      385\n"
     62  "1 5  8      561\n"
     63  "5 3  9      484\n"
     64  "7 5  10     904\n"
     65  "3 6  11     47\n"
     66  "7 6  12     888\n"
     67  "3 0  13     747\n"
     68  "6 1  14     310\n",
     69
     70  "@nodes\n"
     71  "label\n"
     72  "0\n"
     73  "1\n"
     74  "2\n"
     75  "3\n"
     76  "4\n"
     77  "5\n"
     78  "6\n"
     79  "7\n"
     80  "@edges\n"
     81  "     label  weight\n"
     82  "2 5  0      710\n"
     83  "0 5  1      241\n"
     84  "2 4  2      856\n"
     85  "2 6  3      762\n"
     86  "4 1  4      747\n"
     87  "6 1  5      962\n"
     88  "4 7  6      723\n"
     89  "1 7  7      661\n"
     90  "2 3  8      376\n"
     91  "1 0  9      416\n"
     92  "6 7  10     391\n",
     93
     94  "@nodes\n"
     95  "label\n"
     96  "0\n"
     97  "1\n"
     98  "2\n"
     99  "3\n"
     100  "4\n"
     101  "5\n"
     102  "6\n"
     103  "7\n"
     104  "@edges\n"
     105  "     label  weight\n"
     106  "6 2  0      553\n"
     107  "0 7  1      653\n"
     108  "6 3  2      22\n"
     109  "4 7  3      846\n"
     110  "7 2  4      981\n"
     111  "7 6  5      250\n"
     112  "5 2  6      539\n",
     113
     114  "@nodes\n"
     115  "label\n"
     116  "0\n"
     117  "@edges\n"
     118  "     label  weight\n"
     119  "0 0  0      100\n"
     120};
     121
     122void checkMaxFractionalMatchingCompile()
     123{
     124  typedef concepts::Graph Graph;
     125  typedef Graph::Node Node;
     126  typedef Graph::Edge Edge;
     127
     128  Graph g;
     129  Node n;
     130  Edge e;
     131
     132  MaxFractionalMatching<Graph> mat_test(g);
     133  const MaxFractionalMatching<Graph>&
     134    const_mat_test = mat_test;
     135
     136  mat_test.init();
     137  mat_test.start();
     138  mat_test.start(true);
     139  mat_test.startPerfect();
     140  mat_test.startPerfect(true);
     141  mat_test.run();
     142  mat_test.run(true);
     143  mat_test.runPerfect();
     144  mat_test.runPerfect(true);
     145
     146  const_mat_test.matchingSize();
     147  const_mat_test.matching(e);
     148  const_mat_test.matching(n);
     149  const MaxFractionalMatching<Graph>::MatchingMap& mmap =
     150    const_mat_test.matchingMap();
     151  e = mmap[n];
     152
     153  const_mat_test.barrier(n);
     154}
     155
     156void checkMaxWeightedFractionalMatchingCompile()
     157{
     158  typedef concepts::Graph Graph;
     159  typedef Graph::Node Node;
     160  typedef Graph::Edge Edge;
     161  typedef Graph::EdgeMap<int> WeightMap;
     162
     163  Graph g;
     164  Node n;
     165  Edge e;
     166  WeightMap w(g);
     167
     168  MaxWeightedFractionalMatching<Graph> mat_test(g, w);
     169  const MaxWeightedFractionalMatching<Graph>&
     170    const_mat_test = mat_test;
     171
     172  mat_test.init();
     173  mat_test.start();
     174  mat_test.run();
     175
     176  const_mat_test.matchingWeight();
     177  const_mat_test.matchingSize();
     178  const_mat_test.matching(e);
     179  const_mat_test.matching(n);
     180  const MaxWeightedFractionalMatching<Graph>::MatchingMap& mmap =
     181    const_mat_test.matchingMap();
     182  e = mmap[n];
     183
     184  const_mat_test.dualValue();
     185  const_mat_test.nodeValue(n);
     186}
     187
     188void checkMaxWeightedPerfectFractionalMatchingCompile()
     189{
     190  typedef concepts::Graph Graph;
     191  typedef Graph::Node Node;
     192  typedef Graph::Edge Edge;
     193  typedef Graph::EdgeMap<int> WeightMap;
     194
     195  Graph g;
     196  Node n;
     197  Edge e;
     198  WeightMap w(g);
     199
     200  MaxWeightedPerfectFractionalMatching<Graph> mat_test(g, w);
     201  const MaxWeightedPerfectFractionalMatching<Graph>&
     202    const_mat_test = mat_test;
     203
     204  mat_test.init();
     205  mat_test.start();
     206  mat_test.run();
     207
     208  const_mat_test.matchingWeight();
     209  const_mat_test.matching(e);
     210  const_mat_test.matching(n);
     211  const MaxWeightedPerfectFractionalMatching<Graph>::MatchingMap& mmap =
     212    const_mat_test.matchingMap();
     213  e = mmap[n];
     214
     215  const_mat_test.dualValue();
     216  const_mat_test.nodeValue(n);
     217}
     218
     219void checkFractionalMatching(const SmartGraph& graph,
     220                             const MaxFractionalMatching<SmartGraph>& mfm,
     221                             bool allow_loops = true) {
     222  int pv = 0;
     223  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     224    int indeg = 0;
     225    for (InArcIt a(graph, n); a != INVALID; ++a) {
     226      if (mfm.matching(graph.source(a)) == a) {
     227        ++indeg;
     228      }
     229    }
     230    if (mfm.matching(n) != INVALID) {
     231      check(indeg == 1, "Invalid matching");
     232      ++pv;
     233    } else {
     234      check(indeg == 0, "Invalid matching");
     235    }
     236  }
     237  check(pv == mfm.matchingSize(), "Wrong matching size");
     238
     239  SmartGraph::NodeMap<bool> processed(graph, false);
     240  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     241    if (processed[n]) continue;
     242    processed[n] = true;
     243    if (mfm.matching(n) == INVALID) continue;
     244    int num = 1;
     245    Node v = graph.target(mfm.matching(n));
     246    while (v != n) {
     247      processed[v] = true;
     248      ++num;
     249      v = graph.target(mfm.matching(v));
     250    }
     251    check(num == 2 || num % 2 == 1, "Wrong cycle size");
     252    check(allow_loops || num != 1, "Wrong cycle size");
     253  }
     254
     255  int anum = 0, bnum = 0;
     256  SmartGraph::NodeMap<bool> neighbours(graph, false);
     257  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     258    if (!mfm.barrier(n)) continue;
     259    ++anum;
     260    for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) {
     261      Node u = graph.source(a);
     262      if (!allow_loops && u == n) continue;
     263      if (!neighbours[u]) {
     264        neighbours[u] = true;
     265        ++bnum;
     266      }
     267    }
     268  }
     269  check(anum - bnum + mfm.matchingSize() == countNodes(graph),
     270        "Wrong barrier");
     271}
     272
     273void checkPerfectFractionalMatching(const SmartGraph& graph,
     274                             const MaxFractionalMatching<SmartGraph>& mfm,
     275                             bool perfect, bool allow_loops = true) {
     276  if (perfect) {
     277    for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     278      int indeg = 0;
     279      for (InArcIt a(graph, n); a != INVALID; ++a) {
     280        if (mfm.matching(graph.source(a)) == a) {
     281          ++indeg;
     282        }
     283      }
     284      check(mfm.matching(n) != INVALID, "Invalid matching");
     285      check(indeg == 1, "Invalid matching");
     286    }
     287  } else {
     288    int anum = 0, bnum = 0;
     289    SmartGraph::NodeMap<bool> neighbours(graph, false);
     290    for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     291      if (!mfm.barrier(n)) continue;
     292      ++anum;
     293      for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) {
     294        Node u = graph.source(a);
     295        if (!allow_loops && u == n) continue;
     296        if (!neighbours[u]) {
     297          neighbours[u] = true;
     298          ++bnum;
     299        }
     300      }
     301    }
     302    check(anum - bnum > 0, "Wrong barrier");
     303  }
     304}
     305
     306void checkWeightedFractionalMatching(const SmartGraph& graph,
     307                   const SmartGraph::EdgeMap<int>& weight,
     308                   const MaxWeightedFractionalMatching<SmartGraph>& mwfm,
     309                   bool allow_loops = true) {
     310  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
     311    if (graph.u(e) == graph.v(e) && !allow_loops) continue;
     312    int rw = mwfm.nodeValue(graph.u(e)) + mwfm.nodeValue(graph.v(e))
     313      - weight[e] * mwfm.dualScale;
     314
     315    check(rw >= 0, "Negative reduced weight");
     316    check(rw == 0 || !mwfm.matching(e),
     317          "Non-zero reduced weight on matching edge");
     318  }
     319
     320  int pv = 0;
     321  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     322    int indeg = 0;
     323    for (InArcIt a(graph, n); a != INVALID; ++a) {
     324      if (mwfm.matching(graph.source(a)) == a) {
     325        ++indeg;
     326      }
     327    }
     328    check(indeg <= 1, "Invalid matching");
     329    if (mwfm.matching(n) != INVALID) {
     330      check(mwfm.nodeValue(n) >= 0, "Invalid node value");
     331      check(indeg == 1, "Invalid matching");
     332      pv += weight[mwfm.matching(n)];
     333      SmartGraph::Node o = graph.target(mwfm.matching(n));
     334    } else {
     335      check(mwfm.nodeValue(n) == 0, "Invalid matching");
     336      check(indeg == 0, "Invalid matching");
     337    }
     338  }
     339
     340  int dv = 0;
     341  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     342    dv += mwfm.nodeValue(n);
     343  }
     344
     345  check(pv * mwfm.dualScale == dv * 2, "Wrong duality");
     346
     347  SmartGraph::NodeMap<bool> processed(graph, false);
     348  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     349    if (processed[n]) continue;
     350    processed[n] = true;
     351    if (mwfm.matching(n) == INVALID) continue;
     352    int num = 1;
     353    Node v = graph.target(mwfm.matching(n));
     354    while (v != n) {
     355      processed[v] = true;
     356      ++num;
     357      v = graph.target(mwfm.matching(v));
     358    }
     359    check(num == 2 || num % 2 == 1, "Wrong cycle size");
     360    check(allow_loops || num != 1, "Wrong cycle size");
     361  }
     362
     363  return;
     364}
     365
     366void checkWeightedPerfectFractionalMatching(const SmartGraph& graph,
     367                const SmartGraph::EdgeMap<int>& weight,
     368                const MaxWeightedPerfectFractionalMatching<SmartGraph>& mwpfm,
     369                bool allow_loops = true) {
     370  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
     371    if (graph.u(e) == graph.v(e) && !allow_loops) continue;
     372    int rw = mwpfm.nodeValue(graph.u(e)) + mwpfm.nodeValue(graph.v(e))
     373      - weight[e] * mwpfm.dualScale;
     374
     375    check(rw >= 0, "Negative reduced weight");
     376    check(rw == 0 || !mwpfm.matching(e),
     377          "Non-zero reduced weight on matching edge");
     378  }
     379
     380  int pv = 0;
     381  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     382    int indeg = 0;
     383    for (InArcIt a(graph, n); a != INVALID; ++a) {
     384      if (mwpfm.matching(graph.source(a)) == a) {
     385        ++indeg;
     386      }
     387    }
     388    check(mwpfm.matching(n) != INVALID, "Invalid perfect matching");
     389    check(indeg == 1, "Invalid perfect matching");
     390    pv += weight[mwpfm.matching(n)];
     391    SmartGraph::Node o = graph.target(mwpfm.matching(n));
     392  }
     393
     394  int dv = 0;
     395  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     396    dv += mwpfm.nodeValue(n);
     397  }
     398
     399  check(pv * mwpfm.dualScale == dv * 2, "Wrong duality");
     400
     401  SmartGraph::NodeMap<bool> processed(graph, false);
     402  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     403    if (processed[n]) continue;
     404    processed[n] = true;
     405    if (mwpfm.matching(n) == INVALID) continue;
     406    int num = 1;
     407    Node v = graph.target(mwpfm.matching(n));
     408    while (v != n) {
     409      processed[v] = true;
     410      ++num;
     411      v = graph.target(mwpfm.matching(v));
     412    }
     413    check(num == 2 || num % 2 == 1, "Wrong cycle size");
     414    check(allow_loops || num != 1, "Wrong cycle size");
     415  }
     416
     417  return;
     418}
     419
     420
     421int main() {
     422
     423  for (int i = 0; i < lgfn; ++i) {
     424    SmartGraph graph;
     425    SmartGraph::EdgeMap<int> weight(graph);
     426
     427    istringstream lgfs(lgf[i]);
     428    graphReader(graph, lgfs).
     429      edgeMap("weight", weight).run();
     430
     431    bool perfect_with_loops;
     432    {
     433      MaxFractionalMatching<SmartGraph> mfm(graph, true);
     434      mfm.run();
     435      checkFractionalMatching(graph, mfm, true);
     436      perfect_with_loops = mfm.matchingSize() == countNodes(graph);
     437    }
     438
     439    bool perfect_without_loops;
     440    {
     441      MaxFractionalMatching<SmartGraph> mfm(graph, false);
     442      mfm.run();
     443      checkFractionalMatching(graph, mfm, false);
     444      perfect_without_loops = mfm.matchingSize() == countNodes(graph);
     445    }
     446
     447    {
     448      MaxFractionalMatching<SmartGraph> mfm(graph, true);
     449      bool result = mfm.runPerfect();
     450      checkPerfectFractionalMatching(graph, mfm, result, true);
     451      check(result == perfect_with_loops, "Wrong perfect matching");
     452    }
     453
     454    {
     455      MaxFractionalMatching<SmartGraph> mfm(graph, false);
     456      bool result = mfm.runPerfect();
     457      checkPerfectFractionalMatching(graph, mfm, result, false);
     458      check(result == perfect_without_loops, "Wrong perfect matching");
     459    }
     460
     461    {
     462      MaxWeightedFractionalMatching<SmartGraph> mwfm(graph, weight, true);
     463      mwfm.run();
     464      checkWeightedFractionalMatching(graph, weight, mwfm, true);
     465    }
     466
     467    {
     468      MaxWeightedFractionalMatching<SmartGraph> mwfm(graph, weight, false);
     469      mwfm.run();
     470      checkWeightedFractionalMatching(graph, weight, mwfm, false);
     471    }
     472
     473    {
     474      MaxWeightedPerfectFractionalMatching<SmartGraph> mwpfm(graph, weight,
     475                                                             true);
     476      bool perfect = mwpfm.run();
     477      check(perfect == (mwpfm.matchingSize() == countNodes(graph)),
     478            "Perfect matching found");
     479      check(perfect == perfect_with_loops, "Wrong perfect matching");
     480
     481      if (perfect) {
     482        checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, true);
     483      }
     484    }
     485
     486    {
     487      MaxWeightedPerfectFractionalMatching<SmartGraph> mwpfm(graph, weight,
     488                                                             false);
     489      bool perfect = mwpfm.run();
     490      check(perfect == (mwpfm.matchingSize() == countNodes(graph)),
     491            "Perfect matching found");
     492      check(perfect == perfect_without_loops, "Wrong perfect matching");
     493
     494      if (perfect) {
     495        checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, false);
     496      }
     497    }
     498
     499  }
     500
     501  return 0;
     502}