deba@869: /* -*- mode: C++; indent-tabs-mode: nil; -*-
deba@869: *
deba@869: * This file is a part of LEMON, a generic C++ optimization library.
deba@869: *
alpar@877: * Copyright (C) 2003-2010
deba@869: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
deba@869: * (Egervary Research Group on Combinatorial Optimization, EGRES).
deba@869: *
deba@869: * Permission to use, modify and distribute this software is granted
deba@869: * provided that this copyright notice appears in all copies. For
deba@869: * precise terms see the accompanying LICENSE file.
deba@869: *
deba@869: * This software is provided "AS IS" with no warranty of any kind,
deba@869: * express or implied, and with no claim as to its suitability for any
deba@869: * purpose.
deba@869: *
deba@869: */
deba@869:
deba@869: #ifndef LEMON_FRACTIONAL_MATCHING_H
deba@869: #define LEMON_FRACTIONAL_MATCHING_H
deba@869:
deba@869: #include <vector>
deba@869: #include <queue>
deba@869: #include <set>
deba@869: #include <limits>
deba@869:
deba@869: #include <lemon/core.h>
deba@869: #include <lemon/unionfind.h>
deba@869: #include <lemon/bin_heap.h>
deba@869: #include <lemon/maps.h>
deba@869: #include <lemon/assert.h>
deba@869: #include <lemon/elevator.h>
deba@869:
deba@869: ///\ingroup matching
deba@869: ///\file
deba@869: ///\brief Fractional matching algorithms in general graphs.
deba@869:
deba@869: namespace lemon {
deba@869:
deba@869: /// \brief Default traits class of MaxFractionalMatching class.
deba@869: ///
deba@869: /// Default traits class of MaxFractionalMatching class.
deba@869: /// \tparam GR Graph type.
deba@869: template <typename GR>
deba@869: struct MaxFractionalMatchingDefaultTraits {
deba@869:
deba@869: /// \brief The type of the graph the algorithm runs on.
deba@869: typedef GR Graph;
deba@869:
deba@869: /// \brief The type of the map that stores the matching.
deba@869: ///
deba@869: /// The type of the map that stores the matching arcs.
deba@869: /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
deba@869: typedef typename Graph::template NodeMap<typename GR::Arc> MatchingMap;
deba@869:
deba@869: /// \brief Instantiates a MatchingMap.
deba@869: ///
deba@869: /// This function instantiates a \ref MatchingMap.
deba@869: /// \param graph The graph for which we would like to define
deba@869: /// the matching map.
deba@869: static MatchingMap* createMatchingMap(const Graph& graph) {
deba@869: return new MatchingMap(graph);
deba@869: }
deba@869:
deba@869: /// \brief The elevator type used by MaxFractionalMatching algorithm.
deba@869: ///
deba@869: /// The elevator type used by MaxFractionalMatching algorithm.
deba@869: ///
deba@869: /// \sa Elevator
deba@869: /// \sa LinkedElevator
deba@869: typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
deba@869:
deba@869: /// \brief Instantiates an Elevator.
deba@869: ///
deba@869: /// This function instantiates an \ref Elevator.
deba@869: /// \param graph The graph for which we would like to define
deba@869: /// the elevator.
deba@869: /// \param max_level The maximum level of the elevator.
deba@869: static Elevator* createElevator(const Graph& graph, int max_level) {
deba@869: return new Elevator(graph, max_level);
deba@869: }
deba@869: };
deba@869:
deba@869: /// \ingroup matching
deba@869: ///
deba@869: /// \brief Max cardinality fractional matching
deba@869: ///
deba@869: /// This class provides an implementation of fractional matching
deba@869: /// algorithm based on push-relabel principle.
deba@869: ///
deba@869: /// The maximum cardinality fractional matching is a relaxation of the
deba@869: /// maximum cardinality matching problem where the odd set constraints
deba@869: /// are omitted.
deba@869: /// It can be formulated with the following linear program.
deba@869: /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
deba@869: /// \f[x_e \ge 0\quad \forall e\in E\f]
deba@869: /// \f[\max \sum_{e\in E}x_e\f]
deba@869: /// where \f$\delta(X)\f$ is the set of edges incident to a node in
deba@869: /// \f$X\f$. The result can be represented as the union of a
deba@869: /// matching with one value edges and a set of odd length cycles
deba@869: /// with half value edges.
deba@869: ///
deba@869: /// The algorithm calculates an optimal fractional matching and a
deba@869: /// barrier. The number of adjacents of any node set minus the size
deba@869: /// of node set is a lower bound on the uncovered nodes in the
deba@869: /// graph. For maximum matching a barrier is computed which
deba@869: /// maximizes this difference.
deba@869: ///
deba@869: /// The algorithm can be executed with the run() function. After it
deba@869: /// the matching (the primal solution) and the barrier (the dual
deba@869: /// solution) can be obtained using the query functions.
deba@869: ///
deba@869: /// The primal solution is multiplied by
deba@871: /// \ref MaxFractionalMatching::primalScale "2".
deba@869: ///
deba@869: /// \tparam GR The undirected graph type the algorithm runs on.
deba@869: #ifdef DOXYGEN
deba@869: template <typename GR, typename TR>
deba@869: #else
deba@869: template <typename GR,
deba@869: typename TR = MaxFractionalMatchingDefaultTraits<GR> >
deba@869: #endif
deba@869: class MaxFractionalMatching {
deba@869: public:
deba@869:
deba@869: /// \brief The \ref MaxFractionalMatchingDefaultTraits "traits
deba@869: /// class" of the algorithm.
deba@869: typedef TR Traits;
deba@869: /// The type of the graph the algorithm runs on.
deba@869: typedef typename TR::Graph Graph;
deba@869: /// The type of the matching map.
deba@869: typedef typename TR::MatchingMap MatchingMap;
deba@869: /// The type of the elevator.
deba@869: typedef typename TR::Elevator Elevator;
deba@869:
deba@869: /// \brief Scaling factor for primal solution
deba@869: ///
deba@869: /// Scaling factor for primal solution.
deba@869: static const int primalScale = 2;
deba@869:
deba@869: private:
deba@869:
deba@869: const Graph &_graph;
deba@869: int _node_num;
deba@869: bool _allow_loops;
deba@869: int _empty_level;
deba@869:
deba@869: TEMPLATE_GRAPH_TYPEDEFS(Graph);
deba@869:
deba@869: bool _local_matching;
deba@869: MatchingMap *_matching;
deba@869:
deba@869: bool _local_level;
deba@869: Elevator *_level;
deba@869:
deba@869: typedef typename Graph::template NodeMap<int> InDegMap;
deba@869: InDegMap *_indeg;
deba@869:
deba@869: void createStructures() {
deba@869: _node_num = countNodes(_graph);
deba@869:
deba@869: if (!_matching) {
deba@869: _local_matching = true;
deba@869: _matching = Traits::createMatchingMap(_graph);
deba@869: }
deba@869: if (!_level) {
deba@869: _local_level = true;
deba@869: _level = Traits::createElevator(_graph, _node_num);
deba@869: }
deba@869: if (!_indeg) {
deba@869: _indeg = new InDegMap(_graph);
deba@869: }
deba@869: }
deba@869:
deba@869: void destroyStructures() {
deba@869: if (_local_matching) {
deba@869: delete _matching;
deba@869: }
deba@869: if (_local_level) {
deba@869: delete _level;
deba@869: }
deba@869: if (_indeg) {
deba@869: delete _indeg;
deba@869: }
deba@869: }
deba@869:
deba@869: void postprocessing() {
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: if ((*_indeg)[n] != 0) continue;
deba@869: _indeg->set(n, -1);
deba@869: Node u = n;
deba@869: while ((*_matching)[u] != INVALID) {
deba@869: Node v = _graph.target((*_matching)[u]);
deba@869: _indeg->set(v, -1);
deba@869: Arc a = _graph.oppositeArc((*_matching)[u]);
deba@869: u = _graph.target((*_matching)[v]);
deba@869: _indeg->set(u, -1);
deba@869: _matching->set(v, a);
deba@869: }
deba@869: }
deba@869:
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: if ((*_indeg)[n] != 1) continue;
deba@869: _indeg->set(n, -1);
deba@869:
deba@869: int num = 1;
deba@869: Node u = _graph.target((*_matching)[n]);
deba@869: while (u != n) {
deba@869: _indeg->set(u, -1);
deba@869: u = _graph.target((*_matching)[u]);
deba@869: ++num;
deba@869: }
deba@869: if (num % 2 == 0 && num > 2) {
deba@869: Arc prev = _graph.oppositeArc((*_matching)[n]);
deba@869: Node v = _graph.target((*_matching)[n]);
deba@869: u = _graph.target((*_matching)[v]);
deba@869: _matching->set(v, prev);
deba@869: while (u != n) {
deba@869: prev = _graph.oppositeArc((*_matching)[u]);
deba@869: v = _graph.target((*_matching)[u]);
deba@869: u = _graph.target((*_matching)[v]);
deba@869: _matching->set(v, prev);
deba@869: }
deba@869: }
deba@869: }
deba@869: }
deba@869:
deba@869: public:
deba@869:
deba@869: typedef MaxFractionalMatching Create;
deba@869:
deba@869: ///\name Named Template Parameters
deba@869:
deba@869: ///@{
deba@869:
deba@869: template <typename T>
deba@869: struct SetMatchingMapTraits : public Traits {
deba@869: typedef T MatchingMap;
deba@869: static MatchingMap *createMatchingMap(const Graph&) {
deba@869: LEMON_ASSERT(false, "MatchingMap is not initialized");
deba@869: return 0; // ignore warnings
deba@869: }
deba@869: };
deba@869:
deba@869: /// \brief \ref named-templ-param "Named parameter" for setting
deba@869: /// MatchingMap type
deba@869: ///
deba@869: /// \ref named-templ-param "Named parameter" for setting MatchingMap
deba@869: /// type.
deba@869: template <typename T>
deba@869: struct SetMatchingMap
deba@869: : public MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > {
deba@869: typedef MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > Create;
deba@869: };
deba@869:
deba@869: template <typename T>
deba@869: struct SetElevatorTraits : public Traits {
deba@869: typedef T Elevator;
deba@869: static Elevator *createElevator(const Graph&, int) {
deba@869: LEMON_ASSERT(false, "Elevator is not initialized");
deba@869: return 0; // ignore warnings
deba@869: }
deba@869: };
deba@869:
deba@869: /// \brief \ref named-templ-param "Named parameter" for setting
deba@869: /// Elevator type
deba@869: ///
deba@869: /// \ref named-templ-param "Named parameter" for setting Elevator
deba@869: /// type. If this named parameter is used, then an external
deba@869: /// elevator object must be passed to the algorithm using the
deba@869: /// \ref elevator(Elevator&) "elevator()" function before calling
deba@869: /// \ref run() or \ref init().
deba@869: /// \sa SetStandardElevator
deba@869: template <typename T>
deba@869: struct SetElevator
deba@869: : public MaxFractionalMatching<Graph, SetElevatorTraits<T> > {
deba@869: typedef MaxFractionalMatching<Graph, SetElevatorTraits<T> > Create;
deba@869: };
deba@869:
deba@869: template <typename T>
deba@869: struct SetStandardElevatorTraits : public Traits {
deba@869: typedef T Elevator;
deba@869: static Elevator *createElevator(const Graph& graph, int max_level) {
deba@869: return new Elevator(graph, max_level);
deba@869: }
deba@869: };
deba@869:
deba@869: /// \brief \ref named-templ-param "Named parameter" for setting
deba@869: /// Elevator type with automatic allocation
deba@869: ///
deba@869: /// \ref named-templ-param "Named parameter" for setting Elevator
deba@869: /// type with automatic allocation.
deba@869: /// The Elevator should have standard constructor interface to be
deba@869: /// able to automatically created by the algorithm (i.e. the
deba@869: /// graph and the maximum level should be passed to it).
deba@869: /// However an external elevator object could also be passed to the
deba@869: /// algorithm with the \ref elevator(Elevator&) "elevator()" function
deba@869: /// before calling \ref run() or \ref init().
deba@869: /// \sa SetElevator
deba@869: template <typename T>
deba@869: struct SetStandardElevator
deba@869: : public MaxFractionalMatching<Graph, SetStandardElevatorTraits<T> > {
deba@869: typedef MaxFractionalMatching<Graph,
deba@869: SetStandardElevatorTraits<T> > Create;
deba@869: };
deba@869:
deba@869: /// @}
deba@869:
deba@869: protected:
deba@869:
deba@869: MaxFractionalMatching() {}
deba@869:
deba@869: public:
deba@869:
deba@869: /// \brief Constructor
deba@869: ///
deba@869: /// Constructor.
deba@869: ///
deba@869: MaxFractionalMatching(const Graph &graph, bool allow_loops = true)
deba@869: : _graph(graph), _allow_loops(allow_loops),
deba@869: _local_matching(false), _matching(0),
deba@869: _local_level(false), _level(0), _indeg(0)
deba@869: {}
deba@869:
deba@869: ~MaxFractionalMatching() {
deba@869: destroyStructures();
deba@869: }
deba@869:
deba@869: /// \brief Sets the matching map.
deba@869: ///
deba@869: /// Sets the matching map.
deba@869: /// If you don't use this function before calling \ref run() or
deba@869: /// \ref init(), an instance will be allocated automatically.
deba@869: /// The destructor deallocates this automatically allocated map,
deba@869: /// of course.
deba@869: /// \return <tt>(*this)</tt>
deba@869: MaxFractionalMatching& matchingMap(MatchingMap& map) {
deba@869: if (_local_matching) {
deba@869: delete _matching;
deba@869: _local_matching = false;
deba@869: }
deba@869: _matching = ↦
deba@869: return *this;
deba@869: }
deba@869:
deba@869: /// \brief Sets the elevator used by algorithm.
deba@869: ///
deba@869: /// Sets the elevator used by algorithm.
deba@869: /// If you don't use this function before calling \ref run() or
deba@869: /// \ref init(), an instance will be allocated automatically.
deba@869: /// The destructor deallocates this automatically allocated elevator,
deba@869: /// of course.
deba@869: /// \return <tt>(*this)</tt>
deba@869: MaxFractionalMatching& elevator(Elevator& elevator) {
deba@869: if (_local_level) {
deba@869: delete _level;
deba@869: _local_level = false;
deba@869: }
deba@869: _level = &elevator;
deba@869: return *this;
deba@869: }
deba@869:
deba@869: /// \brief Returns a const reference to the elevator.
deba@869: ///
deba@869: /// Returns a const reference to the elevator.
deba@869: ///
deba@869: /// \pre Either \ref run() or \ref init() must be called before
deba@869: /// using this function.
deba@869: const Elevator& elevator() const {
deba@869: return *_level;
deba@869: }
deba@869:
deba@869: /// \name Execution control
deba@869: /// The simplest way to execute the algorithm is to use one of the
deba@869: /// member functions called \c run(). \n
deba@869: /// If you need more control on the execution, first
deba@869: /// you must call \ref init() and then one variant of the start()
deba@869: /// member.
deba@869:
deba@869: /// @{
deba@869:
deba@869: /// \brief Initializes the internal data structures.
deba@869: ///
deba@869: /// Initializes the internal data structures and sets the initial
deba@869: /// matching.
deba@869: void init() {
deba@869: createStructures();
deba@869:
deba@869: _level->initStart();
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: _indeg->set(n, 0);
deba@869: _matching->set(n, INVALID);
deba@869: _level->initAddItem(n);
deba@869: }
deba@869: _level->initFinish();
deba@869:
deba@869: _empty_level = _node_num;
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: for (OutArcIt a(_graph, n); a != INVALID; ++a) {
deba@869: if (_graph.target(a) == n && !_allow_loops) continue;
deba@869: _matching->set(n, a);
deba@869: Node v = _graph.target((*_matching)[n]);
deba@869: _indeg->set(v, (*_indeg)[v] + 1);
deba@869: break;
deba@869: }
deba@869: }
deba@869:
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: if ((*_indeg)[n] == 0) {
deba@869: _level->activate(n);
deba@869: }
deba@869: }
deba@869: }
deba@869:
deba@869: /// \brief Starts the algorithm and computes a fractional matching
deba@869: ///
deba@869: /// The algorithm computes a maximum fractional matching.
deba@869: ///
deba@869: /// \param postprocess The algorithm computes first a matching
deba@869: /// which is a union of a matching with one value edges, cycles
deba@869: /// with half value edges and even length paths with half value
deba@869: /// edges. If the parameter is true, then after the push-relabel
deba@869: /// algorithm it postprocesses the matching to contain only
deba@869: /// matching edges and half value odd cycles.
deba@869: void start(bool postprocess = true) {
deba@869: Node n;
deba@869: while ((n = _level->highestActive()) != INVALID) {
deba@869: int level = _level->highestActiveLevel();
deba@869: int new_level = _level->maxLevel();
deba@869: for (InArcIt a(_graph, n); a != INVALID; ++a) {
deba@869: Node u = _graph.source(a);
deba@869: if (n == u && !_allow_loops) continue;
deba@869: Node v = _graph.target((*_matching)[u]);
deba@869: if ((*_level)[v] < level) {
deba@869: _indeg->set(v, (*_indeg)[v] - 1);
deba@869: if ((*_indeg)[v] == 0) {
deba@869: _level->activate(v);
deba@869: }
deba@869: _matching->set(u, a);
deba@869: _indeg->set(n, (*_indeg)[n] + 1);
deba@869: _level->deactivate(n);
deba@869: goto no_more_push;
deba@869: } else if (new_level > (*_level)[v]) {
deba@869: new_level = (*_level)[v];
deba@869: }
deba@869: }
deba@869:
deba@869: if (new_level + 1 < _level->maxLevel()) {
deba@869: _level->liftHighestActive(new_level + 1);
deba@869: } else {
deba@869: _level->liftHighestActiveToTop();
deba@869: }
deba@869: if (_level->emptyLevel(level)) {
deba@869: _level->liftToTop(level);
deba@869: }
deba@869: no_more_push:
deba@869: ;
deba@869: }
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: if ((*_matching)[n] == INVALID) continue;
deba@869: Node u = _graph.target((*_matching)[n]);
deba@869: if ((*_indeg)[u] > 1) {
deba@869: _indeg->set(u, (*_indeg)[u] - 1);
deba@869: _matching->set(n, INVALID);
deba@869: }
deba@869: }
deba@869: if (postprocess) {
deba@869: postprocessing();
deba@869: }
deba@869: }
deba@869:
deba@869: /// \brief Starts the algorithm and computes a perfect fractional
deba@869: /// matching
deba@869: ///
deba@869: /// The algorithm computes a perfect fractional matching. If it
deba@869: /// does not exists, then the algorithm returns false and the
deba@869: /// matching is undefined and the barrier.
deba@869: ///
deba@869: /// \param postprocess The algorithm computes first a matching
deba@869: /// which is a union of a matching with one value edges, cycles
deba@869: /// with half value edges and even length paths with half value
deba@869: /// edges. If the parameter is true, then after the push-relabel
deba@869: /// algorithm it postprocesses the matching to contain only
deba@869: /// matching edges and half value odd cycles.
deba@869: bool startPerfect(bool postprocess = true) {
deba@869: Node n;
deba@869: while ((n = _level->highestActive()) != INVALID) {
deba@869: int level = _level->highestActiveLevel();
deba@869: int new_level = _level->maxLevel();
deba@869: for (InArcIt a(_graph, n); a != INVALID; ++a) {
deba@869: Node u = _graph.source(a);
deba@869: if (n == u && !_allow_loops) continue;
deba@869: Node v = _graph.target((*_matching)[u]);
deba@869: if ((*_level)[v] < level) {
deba@869: _indeg->set(v, (*_indeg)[v] - 1);
deba@869: if ((*_indeg)[v] == 0) {
deba@869: _level->activate(v);
deba@869: }
deba@869: _matching->set(u, a);
deba@869: _indeg->set(n, (*_indeg)[n] + 1);
deba@869: _level->deactivate(n);
deba@869: goto no_more_push;
deba@869: } else if (new_level > (*_level)[v]) {
deba@869: new_level = (*_level)[v];
deba@869: }
deba@869: }
deba@869:
deba@869: if (new_level + 1 < _level->maxLevel()) {
deba@869: _level->liftHighestActive(new_level + 1);
deba@869: } else {
deba@869: _level->liftHighestActiveToTop();
deba@869: _empty_level = _level->maxLevel() - 1;
deba@869: return false;
deba@869: }
deba@869: if (_level->emptyLevel(level)) {
deba@869: _level->liftToTop(level);
deba@869: _empty_level = level;
deba@869: return false;
deba@869: }
deba@869: no_more_push:
deba@869: ;
deba@869: }
deba@869: if (postprocess) {
deba@869: postprocessing();
deba@869: }
deba@869: return true;
deba@869: }
deba@869:
deba@869: /// \brief Runs the algorithm
deba@869: ///
deba@869: /// Just a shortcut for the next code:
deba@869: ///\code
deba@869: /// init();
deba@869: /// start();
deba@869: ///\endcode
deba@869: void run(bool postprocess = true) {
deba@869: init();
deba@869: start(postprocess);
deba@869: }
deba@869:
deba@869: /// \brief Runs the algorithm to find a perfect fractional matching
deba@869: ///
deba@869: /// Just a shortcut for the next code:
deba@869: ///\code
deba@869: /// init();
deba@869: /// startPerfect();
deba@869: ///\endcode
deba@869: bool runPerfect(bool postprocess = true) {
deba@869: init();
deba@869: return startPerfect(postprocess);
deba@869: }
deba@869:
deba@869: ///@}
deba@869:
deba@869: /// \name Query Functions
deba@869: /// The result of the %Matching algorithm can be obtained using these
deba@869: /// functions.\n
deba@869: /// Before the use of these functions,
deba@869: /// either run() or start() must be called.
deba@869: ///@{
deba@869:
deba@869:
deba@869: /// \brief Return the number of covered nodes in the matching.
deba@869: ///
deba@869: /// This function returns the number of covered nodes in the matching.
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@869: int matchingSize() const {
deba@869: int num = 0;
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: if ((*_matching)[n] != INVALID) {
deba@869: ++num;
deba@869: }
deba@869: }
deba@869: return num;
deba@869: }
deba@869:
deba@869: /// \brief Returns a const reference to the matching map.
deba@869: ///
deba@869: /// Returns a const reference to the node map storing the found
deba@869: /// fractional matching. This method can be called after
deba@869: /// running the algorithm.
deba@869: ///
deba@869: /// \pre Either \ref run() or \ref init() must be called before
deba@869: /// using this function.
deba@869: const MatchingMap& matchingMap() const {
deba@869: return *_matching;
deba@869: }
deba@869:
deba@869: /// \brief Return \c true if the given edge is in the matching.
deba@869: ///
deba@869: /// This function returns \c true if the given edge is in the
deba@869: /// found matching. The result is scaled by \ref primalScale
deba@869: /// "primal scale".
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@869: int matching(const Edge& edge) const {
deba@869: return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0) +
deba@869: (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
deba@869: }
deba@869:
deba@869: /// \brief Return the fractional matching arc (or edge) incident
deba@869: /// to the given node.
deba@869: ///
deba@869: /// This function returns one of the fractional matching arc (or
deba@869: /// edge) incident to the given node in the found matching or \c
deba@869: /// INVALID if the node is not covered by the matching or if the
deba@869: /// node is on an odd length cycle then it is the successor edge
deba@869: /// on the cycle.
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@869: Arc matching(const Node& node) const {
deba@869: return (*_matching)[node];
deba@869: }
deba@869:
deba@869: /// \brief Returns true if the node is in the barrier
deba@869: ///
deba@869: /// The barrier is a subset of the nodes. If the nodes in the
deba@869: /// barrier have less adjacent nodes than the size of the barrier,
deba@869: /// then at least as much nodes cannot be covered as the
deba@869: /// difference of the two subsets.
deba@869: bool barrier(const Node& node) const {
deba@869: return (*_level)[node] >= _empty_level;
deba@869: }
deba@869:
deba@869: /// @}
deba@869:
deba@869: };
deba@869:
deba@869: /// \ingroup matching
deba@869: ///
deba@869: /// \brief Weighted fractional matching in general graphs
deba@869: ///
deba@869: /// This class provides an efficient implementation of fractional
deba@871: /// matching algorithm. The implementation uses priority queues and
deba@871: /// provides \f$O(nm\log n)\f$ time complexity.
deba@869: ///
deba@869: /// The maximum weighted fractional matching is a relaxation of the
deba@869: /// maximum weighted matching problem where the odd set constraints
deba@869: /// are omitted.
deba@869: /// It can be formulated with the following linear program.
deba@869: /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
deba@869: /// \f[x_e \ge 0\quad \forall e\in E\f]
deba@869: /// \f[\max \sum_{e\in E}x_ew_e\f]
deba@869: /// where \f$\delta(X)\f$ is the set of edges incident to a node in
deba@869: /// \f$X\f$. The result must be the union of a matching with one
deba@869: /// value edges and a set of odd length cycles with half value edges.
deba@869: ///
deba@869: /// The algorithm calculates an optimal fractional matching and a
deba@869: /// proof of the optimality. The solution of the dual problem can be
deba@869: /// used to check the result of the algorithm. The dual linear
deba@869: /// problem is the following.
deba@869: /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
deba@869: /// \f[y_u \ge 0 \quad \forall u \in V\f]
deba@871: /// \f[\min \sum_{u \in V}y_u \f]
deba@869: ///
deba@869: /// The algorithm can be executed with the run() function.
deba@869: /// After it the matching (the primal solution) and the dual solution
deba@869: /// can be obtained using the query functions.
deba@869: ///
deba@872: /// The primal solution is multiplied by
deba@872: /// \ref MaxWeightedFractionalMatching::primalScale "2".
deba@872: /// If the value type is integer, then the dual
deba@872: /// solution is scaled by
deba@872: /// \ref MaxWeightedFractionalMatching::dualScale "4".
deba@869: ///
deba@869: /// \tparam GR The undirected graph type the algorithm runs on.
deba@869: /// \tparam WM The type edge weight map. The default type is
deba@869: /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
deba@869: #ifdef DOXYGEN
deba@869: template <typename GR, typename WM>
deba@869: #else
deba@869: template <typename GR,
deba@869: typename WM = typename GR::template EdgeMap<int> >
deba@869: #endif
deba@869: class MaxWeightedFractionalMatching {
deba@869: public:
deba@869:
deba@869: /// The graph type of the algorithm
deba@869: typedef GR Graph;
deba@869: /// The type of the edge weight map
deba@869: typedef WM WeightMap;
deba@869: /// The value type of the edge weights
deba@869: typedef typename WeightMap::Value Value;
deba@869:
deba@869: /// The type of the matching map
deba@869: typedef typename Graph::template NodeMap<typename Graph::Arc>
deba@869: MatchingMap;
deba@869:
deba@869: /// \brief Scaling factor for primal solution
deba@869: ///
deba@872: /// Scaling factor for primal solution.
deba@872: static const int primalScale = 2;
deba@869:
deba@869: /// \brief Scaling factor for dual solution
deba@869: ///
deba@869: /// Scaling factor for dual solution. It is equal to 4 or 1
deba@869: /// according to the value type.
deba@869: static const int dualScale =
deba@869: std::numeric_limits<Value>::is_integer ? 4 : 1;
deba@869:
deba@869: private:
deba@869:
deba@869: TEMPLATE_GRAPH_TYPEDEFS(Graph);
deba@869:
deba@869: typedef typename Graph::template NodeMap<Value> NodePotential;
deba@869:
deba@869: const Graph& _graph;
deba@869: const WeightMap& _weight;
deba@869:
deba@869: MatchingMap* _matching;
deba@869: NodePotential* _node_potential;
deba@869:
deba@869: int _node_num;
deba@869: bool _allow_loops;
deba@869:
deba@869: enum Status {
deba@869: EVEN = -1, MATCHED = 0, ODD = 1
deba@869: };
deba@869:
deba@869: typedef typename Graph::template NodeMap<Status> StatusMap;
deba@869: StatusMap* _status;
deba@869:
deba@869: typedef typename Graph::template NodeMap<Arc> PredMap;
deba@869: PredMap* _pred;
deba@869:
deba@869: typedef ExtendFindEnum<IntNodeMap> TreeSet;
deba@869:
deba@869: IntNodeMap *_tree_set_index;
deba@869: TreeSet *_tree_set;
deba@869:
deba@869: IntNodeMap *_delta1_index;
deba@869: BinHeap<Value, IntNodeMap> *_delta1;
deba@869:
deba@869: IntNodeMap *_delta2_index;
deba@869: BinHeap<Value, IntNodeMap> *_delta2;
deba@869:
deba@869: IntEdgeMap *_delta3_index;
deba@869: BinHeap<Value, IntEdgeMap> *_delta3;
deba@869:
deba@869: Value _delta_sum;
deba@869:
deba@869: void createStructures() {
deba@869: _node_num = countNodes(_graph);
deba@869:
deba@869: if (!_matching) {
deba@869: _matching = new MatchingMap(_graph);
deba@869: }
deba@869: if (!_node_potential) {
deba@869: _node_potential = new NodePotential(_graph);
deba@869: }
deba@869: if (!_status) {
deba@869: _status = new StatusMap(_graph);
deba@869: }
deba@869: if (!_pred) {
deba@869: _pred = new PredMap(_graph);
deba@869: }
deba@869: if (!_tree_set) {
deba@869: _tree_set_index = new IntNodeMap(_graph);
deba@869: _tree_set = new TreeSet(*_tree_set_index);
deba@869: }
deba@869: if (!_delta1) {
deba@869: _delta1_index = new IntNodeMap(_graph);
deba@869: _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
deba@869: }
deba@869: if (!_delta2) {
deba@869: _delta2_index = new IntNodeMap(_graph);
deba@869: _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
deba@869: }
deba@869: if (!_delta3) {
deba@869: _delta3_index = new IntEdgeMap(_graph);
deba@869: _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
deba@869: }
deba@869: }
deba@869:
deba@869: void destroyStructures() {
deba@869: if (_matching) {
deba@869: delete _matching;
deba@869: }
deba@869: if (_node_potential) {
deba@869: delete _node_potential;
deba@869: }
deba@869: if (_status) {
deba@869: delete _status;
deba@869: }
deba@869: if (_pred) {
deba@869: delete _pred;
deba@869: }
deba@869: if (_tree_set) {
deba@869: delete _tree_set_index;
deba@869: delete _tree_set;
deba@869: }
deba@869: if (_delta1) {
deba@869: delete _delta1_index;
deba@869: delete _delta1;
deba@869: }
deba@869: if (_delta2) {
deba@869: delete _delta2_index;
deba@869: delete _delta2;
deba@869: }
deba@869: if (_delta3) {
deba@869: delete _delta3_index;
deba@869: delete _delta3;
deba@869: }
deba@869: }
deba@869:
deba@869: void matchedToEven(Node node, int tree) {
deba@869: _tree_set->insert(node, tree);
deba@869: _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
deba@869: _delta1->push(node, (*_node_potential)[node]);
deba@869:
deba@869: if (_delta2->state(node) == _delta2->IN_HEAP) {
deba@869: _delta2->erase(node);
deba@869: }
deba@869:
deba@869: for (InArcIt a(_graph, node); a != INVALID; ++a) {
deba@869: Node v = _graph.source(a);
deba@869: Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
deba@869: dualScale * _weight[a];
deba@869: if (node == v) {
deba@869: if (_allow_loops && _graph.direction(a)) {
deba@869: _delta3->push(a, rw / 2);
deba@869: }
deba@869: } else if ((*_status)[v] == EVEN) {
deba@869: _delta3->push(a, rw / 2);
deba@869: } else if ((*_status)[v] == MATCHED) {
deba@869: if (_delta2->state(v) != _delta2->IN_HEAP) {
deba@869: _pred->set(v, a);
deba@869: _delta2->push(v, rw);
deba@869: } else if ((*_delta2)[v] > rw) {
deba@869: _pred->set(v, a);
deba@869: _delta2->decrease(v, rw);
deba@869: }
deba@869: }
deba@869: }
deba@869: }
deba@869:
deba@869: void matchedToOdd(Node node, int tree) {
deba@869: _tree_set->insert(node, tree);
deba@869: _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
deba@869:
deba@869: if (_delta2->state(node) == _delta2->IN_HEAP) {
deba@869: _delta2->erase(node);
deba@869: }
deba@869: }
deba@869:
deba@869: void evenToMatched(Node node, int tree) {
deba@869: _delta1->erase(node);
deba@869: _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
deba@869: Arc min = INVALID;
deba@869: Value minrw = std::numeric_limits<Value>::max();
deba@869: for (InArcIt a(_graph, node); a != INVALID; ++a) {
deba@869: Node v = _graph.source(a);
deba@869: Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
deba@869: dualScale * _weight[a];
deba@869:
deba@869: if (node == v) {
deba@869: if (_allow_loops && _graph.direction(a)) {
deba@869: _delta3->erase(a);
deba@869: }
deba@869: } else if ((*_status)[v] == EVEN) {
deba@869: _delta3->erase(a);
deba@869: if (minrw > rw) {
deba@869: min = _graph.oppositeArc(a);
deba@869: minrw = rw;
deba@869: }
deba@869: } else if ((*_status)[v] == MATCHED) {
deba@869: if ((*_pred)[v] == a) {
deba@869: Arc mina = INVALID;
deba@869: Value minrwa = std::numeric_limits<Value>::max();
deba@869: for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
deba@869: Node va = _graph.target(aa);
deba@869: if ((*_status)[va] != EVEN ||
deba@869: _tree_set->find(va) == tree) continue;
deba@869: Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
deba@869: dualScale * _weight[aa];
deba@869: if (minrwa > rwa) {
deba@869: minrwa = rwa;
deba@869: mina = aa;
deba@869: }
deba@869: }
deba@869: if (mina != INVALID) {
deba@869: _pred->set(v, mina);
deba@869: _delta2->increase(v, minrwa);
deba@869: } else {
deba@869: _pred->set(v, INVALID);
deba@869: _delta2->erase(v);
deba@869: }
deba@869: }
deba@869: }
deba@869: }
deba@869: if (min != INVALID) {
deba@869: _pred->set(node, min);
deba@869: _delta2->push(node, minrw);
deba@869: } else {
deba@869: _pred->set(node, INVALID);
deba@869: }
deba@869: }
deba@869:
deba@869: void oddToMatched(Node node) {
deba@869: _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
deba@869: Arc min = INVALID;
deba@869: Value minrw = std::numeric_limits<Value>::max();
deba@869: for (InArcIt a(_graph, node); a != INVALID; ++a) {
deba@869: Node v = _graph.source(a);
deba@869: if ((*_status)[v] != EVEN) continue;
deba@869: Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
deba@869: dualScale * _weight[a];
deba@869:
deba@869: if (minrw > rw) {
deba@869: min = _graph.oppositeArc(a);
deba@869: minrw = rw;
deba@869: }
deba@869: }
deba@869: if (min != INVALID) {
deba@869: _pred->set(node, min);
deba@869: _delta2->push(node, minrw);
deba@869: } else {
deba@869: _pred->set(node, INVALID);
deba@869: }
deba@869: }
deba@869:
deba@869: void alternatePath(Node even, int tree) {
deba@869: Node odd;
deba@869:
deba@869: _status->set(even, MATCHED);
deba@869: evenToMatched(even, tree);
deba@869:
deba@869: Arc prev = (*_matching)[even];
deba@869: while (prev != INVALID) {
deba@869: odd = _graph.target(prev);
deba@869: even = _graph.target((*_pred)[odd]);
deba@869: _matching->set(odd, (*_pred)[odd]);
deba@869: _status->set(odd, MATCHED);
deba@869: oddToMatched(odd);
deba@869:
deba@869: prev = (*_matching)[even];
deba@869: _status->set(even, MATCHED);
deba@869: _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
deba@869: evenToMatched(even, tree);
deba@869: }
deba@869: }
deba@869:
deba@869: void destroyTree(int tree) {
deba@869: for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
deba@869: if ((*_status)[n] == EVEN) {
deba@869: _status->set(n, MATCHED);
deba@869: evenToMatched(n, tree);
deba@869: } else if ((*_status)[n] == ODD) {
deba@869: _status->set(n, MATCHED);
deba@869: oddToMatched(n);
deba@869: }
deba@869: }
deba@869: _tree_set->eraseClass(tree);
deba@869: }
deba@869:
deba@869:
deba@869: void unmatchNode(const Node& node) {
deba@869: int tree = _tree_set->find(node);
deba@869:
deba@869: alternatePath(node, tree);
deba@869: destroyTree(tree);
deba@869:
deba@869: _matching->set(node, INVALID);
deba@869: }
deba@869:
deba@869:
deba@869: void augmentOnEdge(const Edge& edge) {
deba@869: Node left = _graph.u(edge);
deba@869: int left_tree = _tree_set->find(left);
deba@869:
deba@869: alternatePath(left, left_tree);
deba@869: destroyTree(left_tree);
deba@869: _matching->set(left, _graph.direct(edge, true));
deba@869:
deba@869: Node right = _graph.v(edge);
deba@869: int right_tree = _tree_set->find(right);
deba@869:
deba@869: alternatePath(right, right_tree);
deba@869: destroyTree(right_tree);
deba@869: _matching->set(right, _graph.direct(edge, false));
deba@869: }
deba@869:
deba@869: void augmentOnArc(const Arc& arc) {
deba@869: Node left = _graph.source(arc);
deba@869: _status->set(left, MATCHED);
deba@869: _matching->set(left, arc);
deba@869: _pred->set(left, arc);
deba@869:
deba@869: Node right = _graph.target(arc);
deba@869: int right_tree = _tree_set->find(right);
deba@869:
deba@869: alternatePath(right, right_tree);
deba@869: destroyTree(right_tree);
deba@869: _matching->set(right, _graph.oppositeArc(arc));
deba@869: }
deba@869:
deba@869: void extendOnArc(const Arc& arc) {
deba@869: Node base = _graph.target(arc);
deba@869: int tree = _tree_set->find(base);
deba@869:
deba@869: Node odd = _graph.source(arc);
deba@869: _tree_set->insert(odd, tree);
deba@869: _status->set(odd, ODD);
deba@869: matchedToOdd(odd, tree);
deba@869: _pred->set(odd, arc);
deba@869:
deba@869: Node even = _graph.target((*_matching)[odd]);
deba@869: _tree_set->insert(even, tree);
deba@869: _status->set(even, EVEN);
deba@869: matchedToEven(even, tree);
deba@869: }
deba@869:
deba@869: void cycleOnEdge(const Edge& edge, int tree) {
deba@869: Node nca = INVALID;
deba@869: std::vector<Node> left_path, right_path;
deba@869:
deba@869: {
deba@869: std::set<Node> left_set, right_set;
deba@869: Node left = _graph.u(edge);
deba@869: left_path.push_back(left);
deba@869: left_set.insert(left);
deba@869:
deba@869: Node right = _graph.v(edge);
deba@869: right_path.push_back(right);
deba@869: right_set.insert(right);
deba@869:
deba@869: while (true) {
deba@869:
deba@869: if (left_set.find(right) != left_set.end()) {
deba@869: nca = right;
deba@869: break;
deba@869: }
deba@869:
deba@869: if ((*_matching)[left] == INVALID) break;
deba@869:
deba@869: left = _graph.target((*_matching)[left]);
deba@869: left_path.push_back(left);
deba@869: left = _graph.target((*_pred)[left]);
deba@869: left_path.push_back(left);
deba@869:
deba@869: left_set.insert(left);
deba@869:
deba@869: if (right_set.find(left) != right_set.end()) {
deba@869: nca = left;
deba@869: break;
deba@869: }
deba@869:
deba@869: if ((*_matching)[right] == INVALID) break;
deba@869:
deba@869: right = _graph.target((*_matching)[right]);
deba@869: right_path.push_back(right);
deba@869: right = _graph.target((*_pred)[right]);
deba@869: right_path.push_back(right);
deba@869:
deba@869: right_set.insert(right);
deba@869:
deba@869: }
deba@869:
deba@869: if (nca == INVALID) {
deba@869: if ((*_matching)[left] == INVALID) {
deba@869: nca = right;
deba@869: while (left_set.find(nca) == left_set.end()) {
deba@869: nca = _graph.target((*_matching)[nca]);
deba@869: right_path.push_back(nca);
deba@869: nca = _graph.target((*_pred)[nca]);
deba@869: right_path.push_back(nca);
deba@869: }
deba@869: } else {
deba@869: nca = left;
deba@869: while (right_set.find(nca) == right_set.end()) {
deba@869: nca = _graph.target((*_matching)[nca]);
deba@869: left_path.push_back(nca);
deba@869: nca = _graph.target((*_pred)[nca]);
deba@869: left_path.push_back(nca);
deba@869: }
deba@869: }
deba@869: }
deba@869: }
deba@869:
deba@869: alternatePath(nca, tree);
deba@869: Arc prev;
deba@869:
deba@869: prev = _graph.direct(edge, true);
deba@869: for (int i = 0; left_path[i] != nca; i += 2) {
deba@869: _matching->set(left_path[i], prev);
deba@869: _status->set(left_path[i], MATCHED);
deba@869: evenToMatched(left_path[i], tree);
deba@869:
deba@869: prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
deba@869: _status->set(left_path[i + 1], MATCHED);
deba@869: oddToMatched(left_path[i + 1]);
deba@869: }
deba@869: _matching->set(nca, prev);
deba@869:
deba@869: for (int i = 0; right_path[i] != nca; i += 2) {
deba@869: _status->set(right_path[i], MATCHED);
deba@869: evenToMatched(right_path[i], tree);
deba@869:
deba@869: _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
deba@869: _status->set(right_path[i + 1], MATCHED);
deba@869: oddToMatched(right_path[i + 1]);
deba@869: }
deba@869:
deba@869: destroyTree(tree);
deba@869: }
deba@869:
deba@869: void extractCycle(const Arc &arc) {
deba@869: Node left = _graph.source(arc);
deba@869: Node odd = _graph.target((*_matching)[left]);
deba@869: Arc prev;
deba@869: while (odd != left) {
deba@869: Node even = _graph.target((*_matching)[odd]);
deba@869: prev = (*_matching)[odd];
deba@869: odd = _graph.target((*_matching)[even]);
deba@869: _matching->set(even, _graph.oppositeArc(prev));
deba@869: }
deba@869: _matching->set(left, arc);
deba@869:
deba@869: Node right = _graph.target(arc);
deba@869: int right_tree = _tree_set->find(right);
deba@869: alternatePath(right, right_tree);
deba@869: destroyTree(right_tree);
deba@869: _matching->set(right, _graph.oppositeArc(arc));
deba@869: }
deba@869:
deba@869: public:
deba@869:
deba@869: /// \brief Constructor
deba@869: ///
deba@869: /// Constructor.
deba@869: MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
deba@869: bool allow_loops = true)
deba@869: : _graph(graph), _weight(weight), _matching(0),
deba@869: _node_potential(0), _node_num(0), _allow_loops(allow_loops),
deba@869: _status(0), _pred(0),
deba@869: _tree_set_index(0), _tree_set(0),
deba@869:
deba@869: _delta1_index(0), _delta1(0),
deba@869: _delta2_index(0), _delta2(0),
deba@869: _delta3_index(0), _delta3(0),
deba@869:
deba@869: _delta_sum() {}
deba@869:
deba@869: ~MaxWeightedFractionalMatching() {
deba@869: destroyStructures();
deba@869: }
deba@869:
deba@869: /// \name Execution Control
deba@869: /// The simplest way to execute the algorithm is to use the
deba@869: /// \ref run() member function.
deba@869:
deba@869: ///@{
deba@869:
deba@869: /// \brief Initialize the algorithm
deba@869: ///
deba@869: /// This function initializes the algorithm.
deba@869: void init() {
deba@869: createStructures();
deba@869:
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: (*_delta1_index)[n] = _delta1->PRE_HEAP;
deba@869: (*_delta2_index)[n] = _delta2->PRE_HEAP;
deba@869: }
deba@869: for (EdgeIt e(_graph); e != INVALID; ++e) {
deba@869: (*_delta3_index)[e] = _delta3->PRE_HEAP;
deba@869: }
deba@869:
deba@876: _delta1->clear();
deba@876: _delta2->clear();
deba@876: _delta3->clear();
deba@876: _tree_set->clear();
deba@876:
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: Value max = 0;
deba@869: for (OutArcIt e(_graph, n); e != INVALID; ++e) {
deba@869: if (_graph.target(e) == n && !_allow_loops) continue;
deba@869: if ((dualScale * _weight[e]) / 2 > max) {
deba@869: max = (dualScale * _weight[e]) / 2;
deba@869: }
deba@869: }
deba@869: _node_potential->set(n, max);
deba@869: _delta1->push(n, max);
deba@869:
deba@869: _tree_set->insert(n);
deba@869:
deba@869: _matching->set(n, INVALID);
deba@869: _status->set(n, EVEN);
deba@869: }
deba@869:
deba@869: for (EdgeIt e(_graph); e != INVALID; ++e) {
deba@869: Node left = _graph.u(e);
deba@869: Node right = _graph.v(e);
deba@869: if (left == right && !_allow_loops) continue;
deba@869: _delta3->push(e, ((*_node_potential)[left] +
deba@869: (*_node_potential)[right] -
deba@869: dualScale * _weight[e]) / 2);
deba@869: }
deba@869: }
deba@869:
deba@869: /// \brief Start the algorithm
deba@869: ///
deba@869: /// This function starts the algorithm.
deba@869: ///
deba@869: /// \pre \ref init() must be called before using this function.
deba@869: void start() {
deba@869: enum OpType {
deba@869: D1, D2, D3
deba@869: };
deba@869:
deba@869: int unmatched = _node_num;
deba@869: while (unmatched > 0) {
deba@869: Value d1 = !_delta1->empty() ?
deba@869: _delta1->prio() : std::numeric_limits<Value>::max();
deba@869:
deba@869: Value d2 = !_delta2->empty() ?
deba@869: _delta2->prio() : std::numeric_limits<Value>::max();
deba@869:
deba@869: Value d3 = !_delta3->empty() ?
deba@869: _delta3->prio() : std::numeric_limits<Value>::max();
deba@869:
deba@869: _delta_sum = d3; OpType ot = D3;
deba@869: if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
deba@869: if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
deba@869:
deba@869: switch (ot) {
deba@869: case D1:
deba@869: {
deba@869: Node n = _delta1->top();
deba@869: unmatchNode(n);
deba@869: --unmatched;
deba@869: }
deba@869: break;
deba@869: case D2:
deba@869: {
deba@869: Node n = _delta2->top();
deba@869: Arc a = (*_pred)[n];
deba@869: if ((*_matching)[n] == INVALID) {
deba@869: augmentOnArc(a);
deba@869: --unmatched;
deba@869: } else {
deba@869: Node v = _graph.target((*_matching)[n]);
deba@869: if ((*_matching)[n] !=
deba@869: _graph.oppositeArc((*_matching)[v])) {
deba@869: extractCycle(a);
deba@869: --unmatched;
deba@869: } else {
deba@869: extendOnArc(a);
deba@869: }
deba@869: }
deba@869: } break;
deba@869: case D3:
deba@869: {
deba@869: Edge e = _delta3->top();
deba@869:
deba@869: Node left = _graph.u(e);
deba@869: Node right = _graph.v(e);
deba@869:
deba@869: int left_tree = _tree_set->find(left);
deba@869: int right_tree = _tree_set->find(right);
deba@869:
deba@869: if (left_tree == right_tree) {
deba@869: cycleOnEdge(e, left_tree);
deba@869: --unmatched;
deba@869: } else {
deba@869: augmentOnEdge(e);
deba@869: unmatched -= 2;
deba@869: }
deba@869: } break;
deba@869: }
deba@869: }
deba@869: }
deba@869:
deba@869: /// \brief Run the algorithm.
deba@869: ///
deba@871: /// This method runs the \c %MaxWeightedFractionalMatching algorithm.
deba@869: ///
deba@869: /// \note mwfm.run() is just a shortcut of the following code.
deba@869: /// \code
deba@869: /// mwfm.init();
deba@869: /// mwfm.start();
deba@869: /// \endcode
deba@869: void run() {
deba@869: init();
deba@869: start();
deba@869: }
deba@869:
deba@869: /// @}
deba@869:
deba@869: /// \name Primal Solution
deba@869: /// Functions to get the primal solution, i.e. the maximum weighted
deba@869: /// matching.\n
deba@869: /// Either \ref run() or \ref start() function should be called before
deba@869: /// using them.
deba@869:
deba@869: /// @{
deba@869:
deba@869: /// \brief Return the weight of the matching.
deba@869: ///
deba@869: /// This function returns the weight of the found matching. This
deba@869: /// value is scaled by \ref primalScale "primal scale".
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@869: Value matchingWeight() const {
deba@869: Value sum = 0;
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: if ((*_matching)[n] != INVALID) {
deba@869: sum += _weight[(*_matching)[n]];
deba@869: }
deba@869: }
deba@869: return sum * primalScale / 2;
deba@869: }
deba@869:
deba@869: /// \brief Return the number of covered nodes in the matching.
deba@869: ///
deba@869: /// This function returns the number of covered nodes in the matching.
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@869: int matchingSize() const {
deba@869: int num = 0;
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: if ((*_matching)[n] != INVALID) {
deba@869: ++num;
deba@869: }
deba@869: }
deba@869: return num;
deba@869: }
deba@869:
deba@869: /// \brief Return \c true if the given edge is in the matching.
deba@869: ///
deba@869: /// This function returns \c true if the given edge is in the
deba@869: /// found matching. The result is scaled by \ref primalScale
deba@869: /// "primal scale".
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@872: int matching(const Edge& edge) const {
deba@872: return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
deba@872: + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
deba@869: }
deba@869:
deba@869: /// \brief Return the fractional matching arc (or edge) incident
deba@869: /// to the given node.
deba@869: ///
deba@869: /// This function returns one of the fractional matching arc (or
deba@869: /// edge) incident to the given node in the found matching or \c
deba@869: /// INVALID if the node is not covered by the matching or if the
deba@869: /// node is on an odd length cycle then it is the successor edge
deba@869: /// on the cycle.
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@869: Arc matching(const Node& node) const {
deba@869: return (*_matching)[node];
deba@869: }
deba@869:
deba@869: /// \brief Return a const reference to the matching map.
deba@869: ///
deba@869: /// This function returns a const reference to a node map that stores
deba@869: /// the matching arc (or edge) incident to each node.
deba@869: const MatchingMap& matchingMap() const {
deba@869: return *_matching;
deba@869: }
deba@869:
deba@869: /// @}
deba@869:
deba@869: /// \name Dual Solution
deba@869: /// Functions to get the dual solution.\n
deba@869: /// Either \ref run() or \ref start() function should be called before
deba@869: /// using them.
deba@869:
deba@869: /// @{
deba@869:
deba@869: /// \brief Return the value of the dual solution.
deba@869: ///
deba@869: /// This function returns the value of the dual solution.
deba@869: /// It should be equal to the primal value scaled by \ref dualScale
deba@869: /// "dual scale".
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@869: Value dualValue() const {
deba@869: Value sum = 0;
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: sum += nodeValue(n);
deba@869: }
deba@869: return sum;
deba@869: }
deba@869:
deba@869: /// \brief Return the dual value (potential) of the given node.
deba@869: ///
deba@869: /// This function returns the dual value (potential) of the given node.
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@869: Value nodeValue(const Node& n) const {
deba@869: return (*_node_potential)[n];
deba@869: }
deba@869:
deba@869: /// @}
deba@869:
deba@869: };
deba@869:
deba@869: /// \ingroup matching
deba@869: ///
deba@869: /// \brief Weighted fractional perfect matching in general graphs
deba@869: ///
deba@869: /// This class provides an efficient implementation of fractional
deba@871: /// matching algorithm. The implementation uses priority queues and
deba@871: /// provides \f$O(nm\log n)\f$ time complexity.
deba@869: ///
deba@869: /// The maximum weighted fractional perfect matching is a relaxation
deba@869: /// of the maximum weighted perfect matching problem where the odd
deba@869: /// set constraints are omitted.
deba@869: /// It can be formulated with the following linear program.
deba@869: /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
deba@869: /// \f[x_e \ge 0\quad \forall e\in E\f]
deba@869: /// \f[\max \sum_{e\in E}x_ew_e\f]
deba@869: /// where \f$\delta(X)\f$ is the set of edges incident to a node in
deba@869: /// \f$X\f$. The result must be the union of a matching with one
deba@869: /// value edges and a set of odd length cycles with half value edges.
deba@869: ///
deba@869: /// The algorithm calculates an optimal fractional matching and a
deba@869: /// proof of the optimality. The solution of the dual problem can be
deba@869: /// used to check the result of the algorithm. The dual linear
deba@869: /// problem is the following.
deba@869: /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
deba@871: /// \f[\min \sum_{u \in V}y_u \f]
deba@869: ///
deba@869: /// The algorithm can be executed with the run() function.
deba@869: /// After it the matching (the primal solution) and the dual solution
deba@869: /// can be obtained using the query functions.
deba@872: ///
deba@872: /// The primal solution is multiplied by
deba@872: /// \ref MaxWeightedPerfectFractionalMatching::primalScale "2".
deba@872: /// If the value type is integer, then the dual
deba@872: /// solution is scaled by
deba@872: /// \ref MaxWeightedPerfectFractionalMatching::dualScale "4".
deba@869: ///
deba@869: /// \tparam GR The undirected graph type the algorithm runs on.
deba@869: /// \tparam WM The type edge weight map. The default type is
deba@869: /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
deba@869: #ifdef DOXYGEN
deba@869: template <typename GR, typename WM>
deba@869: #else
deba@869: template <typename GR,
deba@869: typename WM = typename GR::template EdgeMap<int> >
deba@869: #endif
deba@869: class MaxWeightedPerfectFractionalMatching {
deba@869: public:
deba@869:
deba@869: /// The graph type of the algorithm
deba@869: typedef GR Graph;
deba@869: /// The type of the edge weight map
deba@869: typedef WM WeightMap;
deba@869: /// The value type of the edge weights
deba@869: typedef typename WeightMap::Value Value;
deba@869:
deba@869: /// The type of the matching map
deba@869: typedef typename Graph::template NodeMap<typename Graph::Arc>
deba@869: MatchingMap;
deba@869:
deba@869: /// \brief Scaling factor for primal solution
deba@869: ///
deba@872: /// Scaling factor for primal solution.
deba@872: static const int primalScale = 2;
deba@869:
deba@869: /// \brief Scaling factor for dual solution
deba@869: ///
deba@869: /// Scaling factor for dual solution. It is equal to 4 or 1
deba@869: /// according to the value type.
deba@869: static const int dualScale =
deba@869: std::numeric_limits<Value>::is_integer ? 4 : 1;
deba@869:
deba@869: private:
deba@869:
deba@869: TEMPLATE_GRAPH_TYPEDEFS(Graph);
deba@869:
deba@869: typedef typename Graph::template NodeMap<Value> NodePotential;
deba@869:
deba@869: const Graph& _graph;
deba@869: const WeightMap& _weight;
deba@869:
deba@869: MatchingMap* _matching;
deba@869: NodePotential* _node_potential;
deba@869:
deba@869: int _node_num;
deba@869: bool _allow_loops;
deba@869:
deba@869: enum Status {
deba@869: EVEN = -1, MATCHED = 0, ODD = 1
deba@869: };
deba@869:
deba@869: typedef typename Graph::template NodeMap<Status> StatusMap;
deba@869: StatusMap* _status;
deba@869:
deba@869: typedef typename Graph::template NodeMap<Arc> PredMap;
deba@869: PredMap* _pred;
deba@869:
deba@869: typedef ExtendFindEnum<IntNodeMap> TreeSet;
deba@869:
deba@869: IntNodeMap *_tree_set_index;
deba@869: TreeSet *_tree_set;
deba@869:
deba@869: IntNodeMap *_delta2_index;
deba@869: BinHeap<Value, IntNodeMap> *_delta2;
deba@869:
deba@869: IntEdgeMap *_delta3_index;
deba@869: BinHeap<Value, IntEdgeMap> *_delta3;
deba@869:
deba@869: Value _delta_sum;
deba@869:
deba@869: void createStructures() {
deba@869: _node_num = countNodes(_graph);
deba@869:
deba@869: if (!_matching) {
deba@869: _matching = new MatchingMap(_graph);
deba@869: }
deba@869: if (!_node_potential) {
deba@869: _node_potential = new NodePotential(_graph);
deba@869: }
deba@869: if (!_status) {
deba@869: _status = new StatusMap(_graph);
deba@869: }
deba@869: if (!_pred) {
deba@869: _pred = new PredMap(_graph);
deba@869: }
deba@869: if (!_tree_set) {
deba@869: _tree_set_index = new IntNodeMap(_graph);
deba@869: _tree_set = new TreeSet(*_tree_set_index);
deba@869: }
deba@869: if (!_delta2) {
deba@869: _delta2_index = new IntNodeMap(_graph);
deba@869: _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
deba@869: }
deba@869: if (!_delta3) {
deba@869: _delta3_index = new IntEdgeMap(_graph);
deba@869: _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
deba@869: }
deba@869: }
deba@869:
deba@869: void destroyStructures() {
deba@869: if (_matching) {
deba@869: delete _matching;
deba@869: }
deba@869: if (_node_potential) {
deba@869: delete _node_potential;
deba@869: }
deba@869: if (_status) {
deba@869: delete _status;
deba@869: }
deba@869: if (_pred) {
deba@869: delete _pred;
deba@869: }
deba@869: if (_tree_set) {
deba@869: delete _tree_set_index;
deba@869: delete _tree_set;
deba@869: }
deba@869: if (_delta2) {
deba@869: delete _delta2_index;
deba@869: delete _delta2;
deba@869: }
deba@869: if (_delta3) {
deba@869: delete _delta3_index;
deba@869: delete _delta3;
deba@869: }
deba@869: }
deba@869:
deba@869: void matchedToEven(Node node, int tree) {
deba@869: _tree_set->insert(node, tree);
deba@869: _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
deba@869:
deba@869: if (_delta2->state(node) == _delta2->IN_HEAP) {
deba@869: _delta2->erase(node);
deba@869: }
deba@869:
deba@869: for (InArcIt a(_graph, node); a != INVALID; ++a) {
deba@869: Node v = _graph.source(a);
deba@869: Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
deba@869: dualScale * _weight[a];
deba@869: if (node == v) {
deba@869: if (_allow_loops && _graph.direction(a)) {
deba@869: _delta3->push(a, rw / 2);
deba@869: }
deba@869: } else if ((*_status)[v] == EVEN) {
deba@869: _delta3->push(a, rw / 2);
deba@869: } else if ((*_status)[v] == MATCHED) {
deba@869: if (_delta2->state(v) != _delta2->IN_HEAP) {
deba@869: _pred->set(v, a);
deba@869: _delta2->push(v, rw);
deba@869: } else if ((*_delta2)[v] > rw) {
deba@869: _pred->set(v, a);
deba@869: _delta2->decrease(v, rw);
deba@869: }
deba@869: }
deba@869: }
deba@869: }
deba@869:
deba@869: void matchedToOdd(Node node, int tree) {
deba@869: _tree_set->insert(node, tree);
deba@869: _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
deba@869:
deba@869: if (_delta2->state(node) == _delta2->IN_HEAP) {
deba@869: _delta2->erase(node);
deba@869: }
deba@869: }
deba@869:
deba@869: void evenToMatched(Node node, int tree) {
deba@869: _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
deba@869: Arc min = INVALID;
deba@869: Value minrw = std::numeric_limits<Value>::max();
deba@869: for (InArcIt a(_graph, node); a != INVALID; ++a) {
deba@869: Node v = _graph.source(a);
deba@869: Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
deba@869: dualScale * _weight[a];
deba@869:
deba@869: if (node == v) {
deba@869: if (_allow_loops && _graph.direction(a)) {
deba@869: _delta3->erase(a);
deba@869: }
deba@869: } else if ((*_status)[v] == EVEN) {
deba@869: _delta3->erase(a);
deba@869: if (minrw > rw) {
deba@869: min = _graph.oppositeArc(a);
deba@869: minrw = rw;
deba@869: }
deba@869: } else if ((*_status)[v] == MATCHED) {
deba@869: if ((*_pred)[v] == a) {
deba@869: Arc mina = INVALID;
deba@869: Value minrwa = std::numeric_limits<Value>::max();
deba@869: for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
deba@869: Node va = _graph.target(aa);
deba@869: if ((*_status)[va] != EVEN ||
deba@869: _tree_set->find(va) == tree) continue;
deba@869: Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
deba@869: dualScale * _weight[aa];
deba@869: if (minrwa > rwa) {
deba@869: minrwa = rwa;
deba@869: mina = aa;
deba@869: }
deba@869: }
deba@869: if (mina != INVALID) {
deba@869: _pred->set(v, mina);
deba@869: _delta2->increase(v, minrwa);
deba@869: } else {
deba@869: _pred->set(v, INVALID);
deba@869: _delta2->erase(v);
deba@869: }
deba@869: }
deba@869: }
deba@869: }
deba@869: if (min != INVALID) {
deba@869: _pred->set(node, min);
deba@869: _delta2->push(node, minrw);
deba@869: } else {
deba@869: _pred->set(node, INVALID);
deba@869: }
deba@869: }
deba@869:
deba@869: void oddToMatched(Node node) {
deba@869: _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
deba@869: Arc min = INVALID;
deba@869: Value minrw = std::numeric_limits<Value>::max();
deba@869: for (InArcIt a(_graph, node); a != INVALID; ++a) {
deba@869: Node v = _graph.source(a);
deba@869: if ((*_status)[v] != EVEN) continue;
deba@869: Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
deba@869: dualScale * _weight[a];
deba@869:
deba@869: if (minrw > rw) {
deba@869: min = _graph.oppositeArc(a);
deba@869: minrw = rw;
deba@869: }
deba@869: }
deba@869: if (min != INVALID) {
deba@869: _pred->set(node, min);
deba@869: _delta2->push(node, minrw);
deba@869: } else {
deba@869: _pred->set(node, INVALID);
deba@869: }
deba@869: }
deba@869:
deba@869: void alternatePath(Node even, int tree) {
deba@869: Node odd;
deba@869:
deba@869: _status->set(even, MATCHED);
deba@869: evenToMatched(even, tree);
deba@869:
deba@869: Arc prev = (*_matching)[even];
deba@869: while (prev != INVALID) {
deba@869: odd = _graph.target(prev);
deba@869: even = _graph.target((*_pred)[odd]);
deba@869: _matching->set(odd, (*_pred)[odd]);
deba@869: _status->set(odd, MATCHED);
deba@869: oddToMatched(odd);
deba@869:
deba@869: prev = (*_matching)[even];
deba@869: _status->set(even, MATCHED);
deba@869: _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
deba@869: evenToMatched(even, tree);
deba@869: }
deba@869: }
deba@869:
deba@869: void destroyTree(int tree) {
deba@869: for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
deba@869: if ((*_status)[n] == EVEN) {
deba@869: _status->set(n, MATCHED);
deba@869: evenToMatched(n, tree);
deba@869: } else if ((*_status)[n] == ODD) {
deba@869: _status->set(n, MATCHED);
deba@869: oddToMatched(n);
deba@869: }
deba@869: }
deba@869: _tree_set->eraseClass(tree);
deba@869: }
deba@869:
deba@869: void augmentOnEdge(const Edge& edge) {
deba@869: Node left = _graph.u(edge);
deba@869: int left_tree = _tree_set->find(left);
deba@869:
deba@869: alternatePath(left, left_tree);
deba@869: destroyTree(left_tree);
deba@869: _matching->set(left, _graph.direct(edge, true));
deba@869:
deba@869: Node right = _graph.v(edge);
deba@869: int right_tree = _tree_set->find(right);
deba@869:
deba@869: alternatePath(right, right_tree);
deba@869: destroyTree(right_tree);
deba@869: _matching->set(right, _graph.direct(edge, false));
deba@869: }
deba@869:
deba@869: void augmentOnArc(const Arc& arc) {
deba@869: Node left = _graph.source(arc);
deba@869: _status->set(left, MATCHED);
deba@869: _matching->set(left, arc);
deba@869: _pred->set(left, arc);
deba@869:
deba@869: Node right = _graph.target(arc);
deba@869: int right_tree = _tree_set->find(right);
deba@869:
deba@869: alternatePath(right, right_tree);
deba@869: destroyTree(right_tree);
deba@869: _matching->set(right, _graph.oppositeArc(arc));
deba@869: }
deba@869:
deba@869: void extendOnArc(const Arc& arc) {
deba@869: Node base = _graph.target(arc);
deba@869: int tree = _tree_set->find(base);
deba@869:
deba@869: Node odd = _graph.source(arc);
deba@869: _tree_set->insert(odd, tree);
deba@869: _status->set(odd, ODD);
deba@869: matchedToOdd(odd, tree);
deba@869: _pred->set(odd, arc);
deba@869:
deba@869: Node even = _graph.target((*_matching)[odd]);
deba@869: _tree_set->insert(even, tree);
deba@869: _status->set(even, EVEN);
deba@869: matchedToEven(even, tree);
deba@869: }
deba@869:
deba@869: void cycleOnEdge(const Edge& edge, int tree) {
deba@869: Node nca = INVALID;
deba@869: std::vector<Node> left_path, right_path;
deba@869:
deba@869: {
deba@869: std::set<Node> left_set, right_set;
deba@869: Node left = _graph.u(edge);
deba@869: left_path.push_back(left);
deba@869: left_set.insert(left);
deba@869:
deba@869: Node right = _graph.v(edge);
deba@869: right_path.push_back(right);
deba@869: right_set.insert(right);
deba@869:
deba@869: while (true) {
deba@869:
deba@869: if (left_set.find(right) != left_set.end()) {
deba@869: nca = right;
deba@869: break;
deba@869: }
deba@869:
deba@869: if ((*_matching)[left] == INVALID) break;
deba@869:
deba@869: left = _graph.target((*_matching)[left]);
deba@869: left_path.push_back(left);
deba@869: left = _graph.target((*_pred)[left]);
deba@869: left_path.push_back(left);
deba@869:
deba@869: left_set.insert(left);
deba@869:
deba@869: if (right_set.find(left) != right_set.end()) {
deba@869: nca = left;
deba@869: break;
deba@869: }
deba@869:
deba@869: if ((*_matching)[right] == INVALID) break;
deba@869:
deba@869: right = _graph.target((*_matching)[right]);
deba@869: right_path.push_back(right);
deba@869: right = _graph.target((*_pred)[right]);
deba@869: right_path.push_back(right);
deba@869:
deba@869: right_set.insert(right);
deba@869:
deba@869: }
deba@869:
deba@869: if (nca == INVALID) {
deba@869: if ((*_matching)[left] == INVALID) {
deba@869: nca = right;
deba@869: while (left_set.find(nca) == left_set.end()) {
deba@869: nca = _graph.target((*_matching)[nca]);
deba@869: right_path.push_back(nca);
deba@869: nca = _graph.target((*_pred)[nca]);
deba@869: right_path.push_back(nca);
deba@869: }
deba@869: } else {
deba@869: nca = left;
deba@869: while (right_set.find(nca) == right_set.end()) {
deba@869: nca = _graph.target((*_matching)[nca]);
deba@869: left_path.push_back(nca);
deba@869: nca = _graph.target((*_pred)[nca]);
deba@869: left_path.push_back(nca);
deba@869: }
deba@869: }
deba@869: }
deba@869: }
deba@869:
deba@869: alternatePath(nca, tree);
deba@869: Arc prev;
deba@869:
deba@869: prev = _graph.direct(edge, true);
deba@869: for (int i = 0; left_path[i] != nca; i += 2) {
deba@869: _matching->set(left_path[i], prev);
deba@869: _status->set(left_path[i], MATCHED);
deba@869: evenToMatched(left_path[i], tree);
deba@869:
deba@869: prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
deba@869: _status->set(left_path[i + 1], MATCHED);
deba@869: oddToMatched(left_path[i + 1]);
deba@869: }
deba@869: _matching->set(nca, prev);
deba@869:
deba@869: for (int i = 0; right_path[i] != nca; i += 2) {
deba@869: _status->set(right_path[i], MATCHED);
deba@869: evenToMatched(right_path[i], tree);
deba@869:
deba@869: _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
deba@869: _status->set(right_path[i + 1], MATCHED);
deba@869: oddToMatched(right_path[i + 1]);
deba@869: }
deba@869:
deba@869: destroyTree(tree);
deba@869: }
deba@869:
deba@869: void extractCycle(const Arc &arc) {
deba@869: Node left = _graph.source(arc);
deba@869: Node odd = _graph.target((*_matching)[left]);
deba@869: Arc prev;
deba@869: while (odd != left) {
deba@869: Node even = _graph.target((*_matching)[odd]);
deba@869: prev = (*_matching)[odd];
deba@869: odd = _graph.target((*_matching)[even]);
deba@869: _matching->set(even, _graph.oppositeArc(prev));
deba@869: }
deba@869: _matching->set(left, arc);
deba@869:
deba@869: Node right = _graph.target(arc);
deba@869: int right_tree = _tree_set->find(right);
deba@869: alternatePath(right, right_tree);
deba@869: destroyTree(right_tree);
deba@869: _matching->set(right, _graph.oppositeArc(arc));
deba@869: }
deba@869:
deba@869: public:
deba@869:
deba@869: /// \brief Constructor
deba@869: ///
deba@869: /// Constructor.
deba@869: MaxWeightedPerfectFractionalMatching(const Graph& graph,
deba@869: const WeightMap& weight,
deba@869: bool allow_loops = true)
deba@869: : _graph(graph), _weight(weight), _matching(0),
deba@869: _node_potential(0), _node_num(0), _allow_loops(allow_loops),
deba@869: _status(0), _pred(0),
deba@869: _tree_set_index(0), _tree_set(0),
deba@869:
deba@869: _delta2_index(0), _delta2(0),
deba@869: _delta3_index(0), _delta3(0),
deba@869:
deba@869: _delta_sum() {}
deba@869:
deba@869: ~MaxWeightedPerfectFractionalMatching() {
deba@869: destroyStructures();
deba@869: }
deba@869:
deba@869: /// \name Execution Control
deba@869: /// The simplest way to execute the algorithm is to use the
deba@869: /// \ref run() member function.
deba@869:
deba@869: ///@{
deba@869:
deba@869: /// \brief Initialize the algorithm
deba@869: ///
deba@869: /// This function initializes the algorithm.
deba@869: void init() {
deba@869: createStructures();
deba@869:
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: (*_delta2_index)[n] = _delta2->PRE_HEAP;
deba@869: }
deba@869: for (EdgeIt e(_graph); e != INVALID; ++e) {
deba@869: (*_delta3_index)[e] = _delta3->PRE_HEAP;
deba@869: }
deba@869:
deba@876: _delta2->clear();
deba@876: _delta3->clear();
deba@876: _tree_set->clear();
deba@876:
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: Value max = - std::numeric_limits<Value>::max();
deba@869: for (OutArcIt e(_graph, n); e != INVALID; ++e) {
deba@869: if (_graph.target(e) == n && !_allow_loops) continue;
deba@869: if ((dualScale * _weight[e]) / 2 > max) {
deba@869: max = (dualScale * _weight[e]) / 2;
deba@869: }
deba@869: }
deba@869: _node_potential->set(n, max);
deba@869:
deba@869: _tree_set->insert(n);
deba@869:
deba@869: _matching->set(n, INVALID);
deba@869: _status->set(n, EVEN);
deba@869: }
deba@869:
deba@869: for (EdgeIt e(_graph); e != INVALID; ++e) {
deba@869: Node left = _graph.u(e);
deba@869: Node right = _graph.v(e);
deba@869: if (left == right && !_allow_loops) continue;
deba@869: _delta3->push(e, ((*_node_potential)[left] +
deba@869: (*_node_potential)[right] -
deba@869: dualScale * _weight[e]) / 2);
deba@869: }
deba@869: }
deba@869:
deba@869: /// \brief Start the algorithm
deba@869: ///
deba@869: /// This function starts the algorithm.
deba@869: ///
deba@869: /// \pre \ref init() must be called before using this function.
deba@869: bool start() {
deba@869: enum OpType {
deba@869: D2, D3
deba@869: };
deba@869:
deba@869: int unmatched = _node_num;
deba@869: while (unmatched > 0) {
deba@869: Value d2 = !_delta2->empty() ?
deba@869: _delta2->prio() : std::numeric_limits<Value>::max();
deba@869:
deba@869: Value d3 = !_delta3->empty() ?
deba@869: _delta3->prio() : std::numeric_limits<Value>::max();
deba@869:
deba@869: _delta_sum = d3; OpType ot = D3;
deba@869: if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
deba@869:
deba@869: if (_delta_sum == std::numeric_limits<Value>::max()) {
deba@869: return false;
deba@869: }
deba@869:
deba@869: switch (ot) {
deba@869: case D2:
deba@869: {
deba@869: Node n = _delta2->top();
deba@869: Arc a = (*_pred)[n];
deba@869: if ((*_matching)[n] == INVALID) {
deba@869: augmentOnArc(a);
deba@869: --unmatched;
deba@869: } else {
deba@869: Node v = _graph.target((*_matching)[n]);
deba@869: if ((*_matching)[n] !=
deba@869: _graph.oppositeArc((*_matching)[v])) {
deba@869: extractCycle(a);
deba@869: --unmatched;
deba@869: } else {
deba@869: extendOnArc(a);
deba@869: }
deba@869: }
deba@869: } break;
deba@869: case D3:
deba@869: {
deba@869: Edge e = _delta3->top();
deba@869:
deba@869: Node left = _graph.u(e);
deba@869: Node right = _graph.v(e);
deba@869:
deba@869: int left_tree = _tree_set->find(left);
deba@869: int right_tree = _tree_set->find(right);
deba@869:
deba@869: if (left_tree == right_tree) {
deba@869: cycleOnEdge(e, left_tree);
deba@869: --unmatched;
deba@869: } else {
deba@869: augmentOnEdge(e);
deba@869: unmatched -= 2;
deba@869: }
deba@869: } break;
deba@869: }
deba@869: }
deba@869: return true;
deba@869: }
deba@869:
deba@869: /// \brief Run the algorithm.
deba@869: ///
alpar@877: /// This method runs the \c %MaxWeightedPerfectFractionalMatching
deba@871: /// algorithm.
deba@869: ///
deba@869: /// \note mwfm.run() is just a shortcut of the following code.
deba@869: /// \code
deba@869: /// mwpfm.init();
deba@869: /// mwpfm.start();
deba@869: /// \endcode
deba@869: bool run() {
deba@869: init();
deba@869: return start();
deba@869: }
deba@869:
deba@869: /// @}
deba@869:
deba@869: /// \name Primal Solution
deba@869: /// Functions to get the primal solution, i.e. the maximum weighted
deba@869: /// matching.\n
deba@869: /// Either \ref run() or \ref start() function should be called before
deba@869: /// using them.
deba@869:
deba@869: /// @{
deba@869:
deba@869: /// \brief Return the weight of the matching.
deba@869: ///
deba@869: /// This function returns the weight of the found matching. This
deba@869: /// value is scaled by \ref primalScale "primal scale".
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@869: Value matchingWeight() const {
deba@869: Value sum = 0;
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: if ((*_matching)[n] != INVALID) {
deba@869: sum += _weight[(*_matching)[n]];
deba@869: }
deba@869: }
deba@869: return sum * primalScale / 2;
deba@869: }
deba@869:
deba@869: /// \brief Return the number of covered nodes in the matching.
deba@869: ///
deba@869: /// This function returns the number of covered nodes in the matching.
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@869: int matchingSize() const {
deba@869: int num = 0;
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: if ((*_matching)[n] != INVALID) {
deba@869: ++num;
deba@869: }
deba@869: }
deba@869: return num;
deba@869: }
deba@869:
deba@869: /// \brief Return \c true if the given edge is in the matching.
deba@869: ///
deba@869: /// This function returns \c true if the given edge is in the
deba@869: /// found matching. The result is scaled by \ref primalScale
deba@869: /// "primal scale".
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@872: int matching(const Edge& edge) const {
deba@872: return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
deba@872: + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
deba@869: }
deba@869:
deba@869: /// \brief Return the fractional matching arc (or edge) incident
deba@869: /// to the given node.
deba@869: ///
deba@869: /// This function returns one of the fractional matching arc (or
deba@869: /// edge) incident to the given node in the found matching or \c
deba@869: /// INVALID if the node is not covered by the matching or if the
deba@869: /// node is on an odd length cycle then it is the successor edge
deba@869: /// on the cycle.
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@869: Arc matching(const Node& node) const {
deba@869: return (*_matching)[node];
deba@869: }
deba@869:
deba@869: /// \brief Return a const reference to the matching map.
deba@869: ///
deba@869: /// This function returns a const reference to a node map that stores
deba@869: /// the matching arc (or edge) incident to each node.
deba@869: const MatchingMap& matchingMap() const {
deba@869: return *_matching;
deba@869: }
deba@869:
deba@869: /// @}
deba@869:
deba@869: /// \name Dual Solution
deba@869: /// Functions to get the dual solution.\n
deba@869: /// Either \ref run() or \ref start() function should be called before
deba@869: /// using them.
deba@869:
deba@869: /// @{
deba@869:
deba@869: /// \brief Return the value of the dual solution.
deba@869: ///
deba@869: /// This function returns the value of the dual solution.
deba@869: /// It should be equal to the primal value scaled by \ref dualScale
deba@869: /// "dual scale".
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@869: Value dualValue() const {
deba@869: Value sum = 0;
deba@869: for (NodeIt n(_graph); n != INVALID; ++n) {
deba@869: sum += nodeValue(n);
deba@869: }
deba@869: return sum;
deba@869: }
deba@869:
deba@869: /// \brief Return the dual value (potential) of the given node.
deba@869: ///
deba@869: /// This function returns the dual value (potential) of the given node.
deba@869: ///
deba@869: /// \pre Either run() or start() must be called before using this function.
deba@869: Value nodeValue(const Node& n) const {
deba@869: return (*_node_potential)[n];
deba@869: }
deba@869:
deba@869: /// @}
deba@869:
deba@869: };
deba@869:
deba@869: } //END OF NAMESPACE LEMON
deba@869:
deba@869: #endif //LEMON_FRACTIONAL_MATCHING_H