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: * deba@869: * Copyright (C) 2003-2009 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 deba@869: #include deba@869: #include deba@869: #include deba@869: deba@869: #include deba@869: #include deba@869: #include deba@869: #include deba@869: #include deba@869: #include 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 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 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 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 deba@869: #else deba@869: template > 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 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 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 deba@869: struct SetMatchingMap deba@869: : public MaxFractionalMatching > { deba@869: typedef MaxFractionalMatching > Create; deba@869: }; deba@869: deba@869: template 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 deba@869: struct SetElevator deba@869: : public MaxFractionalMatching > { deba@869: typedef MaxFractionalMatching > Create; deba@869: }; deba@869: deba@869: template 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 deba@869: struct SetStandardElevator deba@869: : public MaxFractionalMatching > { deba@869: typedef MaxFractionalMatching > 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 (*this) 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 (*this) 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". deba@869: #ifdef DOXYGEN deba@869: template deba@869: #else deba@869: template > 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 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::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 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 StatusMap; deba@869: StatusMap* _status; deba@869: deba@869: typedef typename Graph::template NodeMap PredMap; deba@869: PredMap* _pred; deba@869: deba@869: typedef ExtendFindEnum TreeSet; deba@869: deba@869: IntNodeMap *_tree_set_index; deba@869: TreeSet *_tree_set; deba@869: deba@869: IntNodeMap *_delta1_index; deba@869: BinHeap *_delta1; deba@869: deba@869: IntNodeMap *_delta2_index; deba@869: BinHeap *_delta2; deba@869: deba@869: IntEdgeMap *_delta3_index; deba@869: BinHeap *_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(*_delta1_index); deba@869: } deba@869: if (!_delta2) { deba@869: _delta2_index = new IntNodeMap(_graph); deba@869: _delta2 = new BinHeap(*_delta2_index); deba@869: } deba@869: if (!_delta3) { deba@869: _delta3_index = new IntEdgeMap(_graph); deba@869: _delta3 = new BinHeap(*_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::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::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::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 left_path, right_path; deba@869: deba@869: { deba@869: std::set 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@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::max(); deba@869: deba@869: Value d2 = !_delta2->empty() ? deba@869: _delta2->prio() : std::numeric_limits::max(); deba@869: deba@869: Value d3 = !_delta3->empty() ? deba@869: _delta3->prio() : std::numeric_limits::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". deba@869: #ifdef DOXYGEN deba@869: template deba@869: #else deba@869: template > 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 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::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 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 StatusMap; deba@869: StatusMap* _status; deba@869: deba@869: typedef typename Graph::template NodeMap PredMap; deba@869: PredMap* _pred; deba@869: deba@869: typedef ExtendFindEnum TreeSet; deba@869: deba@869: IntNodeMap *_tree_set_index; deba@869: TreeSet *_tree_set; deba@869: deba@869: IntNodeMap *_delta2_index; deba@869: BinHeap *_delta2; deba@869: deba@869: IntEdgeMap *_delta3_index; deba@869: BinHeap *_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(*_delta2_index); deba@869: } deba@869: if (!_delta3) { deba@869: _delta3_index = new IntEdgeMap(_graph); deba@869: _delta3 = new BinHeap(*_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::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::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::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 left_path, right_path; deba@869: deba@869: { deba@869: std::set 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@869: for (NodeIt n(_graph); n != INVALID; ++n) { deba@869: Value max = - std::numeric_limits::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::max(); deba@869: deba@869: Value d3 = !_delta3->empty() ? deba@869: _delta3->prio() : std::numeric_limits::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::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: /// deba@871: /// 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