# HG changeset patch # User Balazs Dezso # Date 1253908296 -7200 # Node ID 636dadefe1e6736f97988b05d6b4e7f95e4c5b9c # Parent 0513ccfea9670a94312cc2bde2ed8595f08361aa Add fractional matching algorithms (#314) diff -r 0513ccfea967 -r 636dadefe1e6 doc/groups.dox --- a/doc/groups.dox Sun Sep 20 21:38:24 2009 +0200 +++ b/doc/groups.dox Fri Sep 25 21:51:36 2009 +0200 @@ -349,7 +349,7 @@ also provide functions to query the minimum cut, which is the dual problem of maximum flow. -\ref Circulation is a preflow push-relabel algorithm implemented directly +\ref Circulation is a preflow push-relabel algorithm implemented directly for finding feasible circulations, which is a somewhat different problem, but it is strongly related to maximum flow. For more information, see \ref Circulation. @@ -470,6 +470,13 @@ - \ref MaxWeightedPerfectMatching Edmond's blossom shrinking algorithm for calculating maximum weighted perfect matching in general graphs. +- \ref MaxFractionalMatching Push-relabel algorithm for calculating + maximum cardinality fractional matching in general graphs. +- \ref MaxWeightedFractionalMatching Augmenting path algorithm for calculating + maximum weighted fractional matching in general graphs. +- \ref MaxWeightedPerfectFractionalMatching + Augmenting path algorithm for calculating maximum weighted + perfect fractional matching in general graphs. \image html bipartite_matching.png \image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth diff -r 0513ccfea967 -r 636dadefe1e6 lemon/Makefile.am --- a/lemon/Makefile.am Sun Sep 20 21:38:24 2009 +0200 +++ b/lemon/Makefile.am Fri Sep 25 21:51:36 2009 +0200 @@ -81,6 +81,7 @@ lemon/euler.h \ lemon/fib_heap.h \ lemon/fourary_heap.h \ + lemon/fractional_matching.h \ lemon/full_graph.h \ lemon/glpk.h \ lemon/gomory_hu.h \ diff -r 0513ccfea967 -r 636dadefe1e6 lemon/fractional_matching.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lemon/fractional_matching.h Fri Sep 25 21:51:36 2009 +0200 @@ -0,0 +1,2135 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2009 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef LEMON_FRACTIONAL_MATCHING_H +#define LEMON_FRACTIONAL_MATCHING_H + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +///\ingroup matching +///\file +///\brief Fractional matching algorithms in general graphs. + +namespace lemon { + + /// \brief Default traits class of MaxFractionalMatching class. + /// + /// Default traits class of MaxFractionalMatching class. + /// \tparam GR Graph type. + template + struct MaxFractionalMatchingDefaultTraits { + + /// \brief The type of the graph the algorithm runs on. + typedef GR Graph; + + /// \brief The type of the map that stores the matching. + /// + /// The type of the map that stores the matching arcs. + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept. + typedef typename Graph::template NodeMap MatchingMap; + + /// \brief Instantiates a MatchingMap. + /// + /// This function instantiates a \ref MatchingMap. + /// \param graph The graph for which we would like to define + /// the matching map. + static MatchingMap* createMatchingMap(const Graph& graph) { + return new MatchingMap(graph); + } + + /// \brief The elevator type used by MaxFractionalMatching algorithm. + /// + /// The elevator type used by MaxFractionalMatching algorithm. + /// + /// \sa Elevator + /// \sa LinkedElevator + typedef LinkedElevator Elevator; + + /// \brief Instantiates an Elevator. + /// + /// This function instantiates an \ref Elevator. + /// \param graph The graph for which we would like to define + /// the elevator. + /// \param max_level The maximum level of the elevator. + static Elevator* createElevator(const Graph& graph, int max_level) { + return new Elevator(graph, max_level); + } + }; + + /// \ingroup matching + /// + /// \brief Max cardinality fractional matching + /// + /// This class provides an implementation of fractional matching + /// algorithm based on push-relabel principle. + /// + /// The maximum cardinality fractional matching is a relaxation of the + /// maximum cardinality matching problem where the odd set constraints + /// are omitted. + /// It can be formulated with the following linear program. + /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f] + /// \f[x_e \ge 0\quad \forall e\in E\f] + /// \f[\max \sum_{e\in E}x_e\f] + /// where \f$\delta(X)\f$ is the set of edges incident to a node in + /// \f$X\f$. The result can be represented as the union of a + /// matching with one value edges and a set of odd length cycles + /// with half value edges. + /// + /// The algorithm calculates an optimal fractional matching and a + /// barrier. The number of adjacents of any node set minus the size + /// of node set is a lower bound on the uncovered nodes in the + /// graph. For maximum matching a barrier is computed which + /// maximizes this difference. + /// + /// The algorithm can be executed with the run() function. After it + /// the matching (the primal solution) and the barrier (the dual + /// solution) can be obtained using the query functions. + /// + /// The primal solution is multiplied by + /// \ref MaxWeightedMatching::primalScale "2". + /// + /// \tparam GR The undirected graph type the algorithm runs on. +#ifdef DOXYGEN + template +#else + template > +#endif + class MaxFractionalMatching { + public: + + /// \brief The \ref MaxFractionalMatchingDefaultTraits "traits + /// class" of the algorithm. + typedef TR Traits; + /// The type of the graph the algorithm runs on. + typedef typename TR::Graph Graph; + /// The type of the matching map. + typedef typename TR::MatchingMap MatchingMap; + /// The type of the elevator. + typedef typename TR::Elevator Elevator; + + /// \brief Scaling factor for primal solution + /// + /// Scaling factor for primal solution. + static const int primalScale = 2; + + private: + + const Graph &_graph; + int _node_num; + bool _allow_loops; + int _empty_level; + + TEMPLATE_GRAPH_TYPEDEFS(Graph); + + bool _local_matching; + MatchingMap *_matching; + + bool _local_level; + Elevator *_level; + + typedef typename Graph::template NodeMap InDegMap; + InDegMap *_indeg; + + void createStructures() { + _node_num = countNodes(_graph); + + if (!_matching) { + _local_matching = true; + _matching = Traits::createMatchingMap(_graph); + } + if (!_level) { + _local_level = true; + _level = Traits::createElevator(_graph, _node_num); + } + if (!_indeg) { + _indeg = new InDegMap(_graph); + } + } + + void destroyStructures() { + if (_local_matching) { + delete _matching; + } + if (_local_level) { + delete _level; + } + if (_indeg) { + delete _indeg; + } + } + + void postprocessing() { + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_indeg)[n] != 0) continue; + _indeg->set(n, -1); + Node u = n; + while ((*_matching)[u] != INVALID) { + Node v = _graph.target((*_matching)[u]); + _indeg->set(v, -1); + Arc a = _graph.oppositeArc((*_matching)[u]); + u = _graph.target((*_matching)[v]); + _indeg->set(u, -1); + _matching->set(v, a); + } + } + + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_indeg)[n] != 1) continue; + _indeg->set(n, -1); + + int num = 1; + Node u = _graph.target((*_matching)[n]); + while (u != n) { + _indeg->set(u, -1); + u = _graph.target((*_matching)[u]); + ++num; + } + if (num % 2 == 0 && num > 2) { + Arc prev = _graph.oppositeArc((*_matching)[n]); + Node v = _graph.target((*_matching)[n]); + u = _graph.target((*_matching)[v]); + _matching->set(v, prev); + while (u != n) { + prev = _graph.oppositeArc((*_matching)[u]); + v = _graph.target((*_matching)[u]); + u = _graph.target((*_matching)[v]); + _matching->set(v, prev); + } + } + } + } + + public: + + typedef MaxFractionalMatching Create; + + ///\name Named Template Parameters + + ///@{ + + template + struct SetMatchingMapTraits : public Traits { + typedef T MatchingMap; + static MatchingMap *createMatchingMap(const Graph&) { + LEMON_ASSERT(false, "MatchingMap is not initialized"); + return 0; // ignore warnings + } + }; + + /// \brief \ref named-templ-param "Named parameter" for setting + /// MatchingMap type + /// + /// \ref named-templ-param "Named parameter" for setting MatchingMap + /// type. + template + struct SetMatchingMap + : public MaxFractionalMatching > { + typedef MaxFractionalMatching > Create; + }; + + template + struct SetElevatorTraits : public Traits { + typedef T Elevator; + static Elevator *createElevator(const Graph&, int) { + LEMON_ASSERT(false, "Elevator is not initialized"); + return 0; // ignore warnings + } + }; + + /// \brief \ref named-templ-param "Named parameter" for setting + /// Elevator type + /// + /// \ref named-templ-param "Named parameter" for setting Elevator + /// type. If this named parameter is used, then an external + /// elevator object must be passed to the algorithm using the + /// \ref elevator(Elevator&) "elevator()" function before calling + /// \ref run() or \ref init(). + /// \sa SetStandardElevator + template + struct SetElevator + : public MaxFractionalMatching > { + typedef MaxFractionalMatching > Create; + }; + + template + struct SetStandardElevatorTraits : public Traits { + typedef T Elevator; + static Elevator *createElevator(const Graph& graph, int max_level) { + return new Elevator(graph, max_level); + } + }; + + /// \brief \ref named-templ-param "Named parameter" for setting + /// Elevator type with automatic allocation + /// + /// \ref named-templ-param "Named parameter" for setting Elevator + /// type with automatic allocation. + /// The Elevator should have standard constructor interface to be + /// able to automatically created by the algorithm (i.e. the + /// graph and the maximum level should be passed to it). + /// However an external elevator object could also be passed to the + /// algorithm with the \ref elevator(Elevator&) "elevator()" function + /// before calling \ref run() or \ref init(). + /// \sa SetElevator + template + struct SetStandardElevator + : public MaxFractionalMatching > { + typedef MaxFractionalMatching > Create; + }; + + /// @} + + protected: + + MaxFractionalMatching() {} + + public: + + /// \brief Constructor + /// + /// Constructor. + /// + MaxFractionalMatching(const Graph &graph, bool allow_loops = true) + : _graph(graph), _allow_loops(allow_loops), + _local_matching(false), _matching(0), + _local_level(false), _level(0), _indeg(0) + {} + + ~MaxFractionalMatching() { + destroyStructures(); + } + + /// \brief Sets the matching map. + /// + /// Sets the matching map. + /// If you don't use this function before calling \ref run() or + /// \ref init(), an instance will be allocated automatically. + /// The destructor deallocates this automatically allocated map, + /// of course. + /// \return (*this) + MaxFractionalMatching& matchingMap(MatchingMap& map) { + if (_local_matching) { + delete _matching; + _local_matching = false; + } + _matching = ↦ + return *this; + } + + /// \brief Sets the elevator used by algorithm. + /// + /// Sets the elevator used by algorithm. + /// If you don't use this function before calling \ref run() or + /// \ref init(), an instance will be allocated automatically. + /// The destructor deallocates this automatically allocated elevator, + /// of course. + /// \return (*this) + MaxFractionalMatching& elevator(Elevator& elevator) { + if (_local_level) { + delete _level; + _local_level = false; + } + _level = &elevator; + return *this; + } + + /// \brief Returns a const reference to the elevator. + /// + /// Returns a const reference to the elevator. + /// + /// \pre Either \ref run() or \ref init() must be called before + /// using this function. + const Elevator& elevator() const { + return *_level; + } + + /// \name Execution control + /// The simplest way to execute the algorithm is to use one of the + /// member functions called \c run(). \n + /// If you need more control on the execution, first + /// you must call \ref init() and then one variant of the start() + /// member. + + /// @{ + + /// \brief Initializes the internal data structures. + /// + /// Initializes the internal data structures and sets the initial + /// matching. + void init() { + createStructures(); + + _level->initStart(); + for (NodeIt n(_graph); n != INVALID; ++n) { + _indeg->set(n, 0); + _matching->set(n, INVALID); + _level->initAddItem(n); + } + _level->initFinish(); + + _empty_level = _node_num; + for (NodeIt n(_graph); n != INVALID; ++n) { + for (OutArcIt a(_graph, n); a != INVALID; ++a) { + if (_graph.target(a) == n && !_allow_loops) continue; + _matching->set(n, a); + Node v = _graph.target((*_matching)[n]); + _indeg->set(v, (*_indeg)[v] + 1); + break; + } + } + + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_indeg)[n] == 0) { + _level->activate(n); + } + } + } + + /// \brief Starts the algorithm and computes a fractional matching + /// + /// The algorithm computes a maximum fractional matching. + /// + /// \param postprocess The algorithm computes first a matching + /// which is a union of a matching with one value edges, cycles + /// with half value edges and even length paths with half value + /// edges. If the parameter is true, then after the push-relabel + /// algorithm it postprocesses the matching to contain only + /// matching edges and half value odd cycles. + void start(bool postprocess = true) { + Node n; + while ((n = _level->highestActive()) != INVALID) { + int level = _level->highestActiveLevel(); + int new_level = _level->maxLevel(); + for (InArcIt a(_graph, n); a != INVALID; ++a) { + Node u = _graph.source(a); + if (n == u && !_allow_loops) continue; + Node v = _graph.target((*_matching)[u]); + if ((*_level)[v] < level) { + _indeg->set(v, (*_indeg)[v] - 1); + if ((*_indeg)[v] == 0) { + _level->activate(v); + } + _matching->set(u, a); + _indeg->set(n, (*_indeg)[n] + 1); + _level->deactivate(n); + goto no_more_push; + } else if (new_level > (*_level)[v]) { + new_level = (*_level)[v]; + } + } + + if (new_level + 1 < _level->maxLevel()) { + _level->liftHighestActive(new_level + 1); + } else { + _level->liftHighestActiveToTop(); + } + if (_level->emptyLevel(level)) { + _level->liftToTop(level); + } + no_more_push: + ; + } + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_matching)[n] == INVALID) continue; + Node u = _graph.target((*_matching)[n]); + if ((*_indeg)[u] > 1) { + _indeg->set(u, (*_indeg)[u] - 1); + _matching->set(n, INVALID); + } + } + if (postprocess) { + postprocessing(); + } + } + + /// \brief Starts the algorithm and computes a perfect fractional + /// matching + /// + /// The algorithm computes a perfect fractional matching. If it + /// does not exists, then the algorithm returns false and the + /// matching is undefined and the barrier. + /// + /// \param postprocess The algorithm computes first a matching + /// which is a union of a matching with one value edges, cycles + /// with half value edges and even length paths with half value + /// edges. If the parameter is true, then after the push-relabel + /// algorithm it postprocesses the matching to contain only + /// matching edges and half value odd cycles. + bool startPerfect(bool postprocess = true) { + Node n; + while ((n = _level->highestActive()) != INVALID) { + int level = _level->highestActiveLevel(); + int new_level = _level->maxLevel(); + for (InArcIt a(_graph, n); a != INVALID; ++a) { + Node u = _graph.source(a); + if (n == u && !_allow_loops) continue; + Node v = _graph.target((*_matching)[u]); + if ((*_level)[v] < level) { + _indeg->set(v, (*_indeg)[v] - 1); + if ((*_indeg)[v] == 0) { + _level->activate(v); + } + _matching->set(u, a); + _indeg->set(n, (*_indeg)[n] + 1); + _level->deactivate(n); + goto no_more_push; + } else if (new_level > (*_level)[v]) { + new_level = (*_level)[v]; + } + } + + if (new_level + 1 < _level->maxLevel()) { + _level->liftHighestActive(new_level + 1); + } else { + _level->liftHighestActiveToTop(); + _empty_level = _level->maxLevel() - 1; + return false; + } + if (_level->emptyLevel(level)) { + _level->liftToTop(level); + _empty_level = level; + return false; + } + no_more_push: + ; + } + if (postprocess) { + postprocessing(); + } + return true; + } + + /// \brief Runs the algorithm + /// + /// Just a shortcut for the next code: + ///\code + /// init(); + /// start(); + ///\endcode + void run(bool postprocess = true) { + init(); + start(postprocess); + } + + /// \brief Runs the algorithm to find a perfect fractional matching + /// + /// Just a shortcut for the next code: + ///\code + /// init(); + /// startPerfect(); + ///\endcode + bool runPerfect(bool postprocess = true) { + init(); + return startPerfect(postprocess); + } + + ///@} + + /// \name Query Functions + /// The result of the %Matching algorithm can be obtained using these + /// functions.\n + /// Before the use of these functions, + /// either run() or start() must be called. + ///@{ + + + /// \brief Return the number of covered nodes in the matching. + /// + /// This function returns the number of covered nodes in the matching. + /// + /// \pre Either run() or start() must be called before using this function. + int matchingSize() const { + int num = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_matching)[n] != INVALID) { + ++num; + } + } + return num; + } + + /// \brief Returns a const reference to the matching map. + /// + /// Returns a const reference to the node map storing the found + /// fractional matching. This method can be called after + /// running the algorithm. + /// + /// \pre Either \ref run() or \ref init() must be called before + /// using this function. + const MatchingMap& matchingMap() const { + return *_matching; + } + + /// \brief Return \c true if the given edge is in the matching. + /// + /// This function returns \c true if the given edge is in the + /// found matching. The result is scaled by \ref primalScale + /// "primal scale". + /// + /// \pre Either run() or start() must be called before using this function. + int matching(const Edge& edge) const { + return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0) + + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0); + } + + /// \brief Return the fractional matching arc (or edge) incident + /// to the given node. + /// + /// This function returns one of the fractional matching arc (or + /// edge) incident to the given node in the found matching or \c + /// INVALID if the node is not covered by the matching or if the + /// node is on an odd length cycle then it is the successor edge + /// on the cycle. + /// + /// \pre Either run() or start() must be called before using this function. + Arc matching(const Node& node) const { + return (*_matching)[node]; + } + + /// \brief Returns true if the node is in the barrier + /// + /// The barrier is a subset of the nodes. If the nodes in the + /// barrier have less adjacent nodes than the size of the barrier, + /// then at least as much nodes cannot be covered as the + /// difference of the two subsets. + bool barrier(const Node& node) const { + return (*_level)[node] >= _empty_level; + } + + /// @} + + }; + + /// \ingroup matching + /// + /// \brief Weighted fractional matching in general graphs + /// + /// This class provides an efficient implementation of fractional + /// matching algorithm. The implementation is based on extensive use + /// of priority queues and provides \f$O(nm\log n)\f$ time + /// complexity. + /// + /// The maximum weighted fractional matching is a relaxation of the + /// maximum weighted matching problem where the odd set constraints + /// are omitted. + /// It can be formulated with the following linear program. + /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f] + /// \f[x_e \ge 0\quad \forall e\in E\f] + /// \f[\max \sum_{e\in E}x_ew_e\f] + /// where \f$\delta(X)\f$ is the set of edges incident to a node in + /// \f$X\f$. The result must be the union of a matching with one + /// value edges and a set of odd length cycles with half value edges. + /// + /// The algorithm calculates an optimal fractional matching and a + /// proof of the optimality. The solution of the dual problem can be + /// used to check the result of the algorithm. The dual linear + /// problem is the following. + /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f] + /// \f[y_u \ge 0 \quad \forall u \in V\f] + /// \f[\min \sum_{u \in V}y_u \f] /// + /// + /// The algorithm can be executed with the run() function. + /// After it the matching (the primal solution) and the dual solution + /// can be obtained using the query functions. + /// + /// If the value type is integer, then the primal and the dual + /// solutions are multiplied by + /// \ref MaxWeightedMatching::primalScale "2" and + /// \ref MaxWeightedMatching::dualScale "4" respectively. + /// + /// \tparam GR The undirected graph type the algorithm runs on. + /// \tparam WM The type edge weight map. The default type is + /// \ref concepts::Graph::EdgeMap "GR::EdgeMap". +#ifdef DOXYGEN + template +#else + template > +#endif + class MaxWeightedFractionalMatching { + public: + + /// The graph type of the algorithm + typedef GR Graph; + /// The type of the edge weight map + typedef WM WeightMap; + /// The value type of the edge weights + typedef typename WeightMap::Value Value; + + /// The type of the matching map + typedef typename Graph::template NodeMap + MatchingMap; + + /// \brief Scaling factor for primal solution + /// + /// Scaling factor for primal solution. It is equal to 2 or 1 + /// according to the value type. + static const int primalScale = + std::numeric_limits::is_integer ? 2 : 1; + + /// \brief Scaling factor for dual solution + /// + /// Scaling factor for dual solution. It is equal to 4 or 1 + /// according to the value type. + static const int dualScale = + std::numeric_limits::is_integer ? 4 : 1; + + private: + + TEMPLATE_GRAPH_TYPEDEFS(Graph); + + typedef typename Graph::template NodeMap NodePotential; + + const Graph& _graph; + const WeightMap& _weight; + + MatchingMap* _matching; + NodePotential* _node_potential; + + int _node_num; + bool _allow_loops; + + enum Status { + EVEN = -1, MATCHED = 0, ODD = 1 + }; + + typedef typename Graph::template NodeMap StatusMap; + StatusMap* _status; + + typedef typename Graph::template NodeMap PredMap; + PredMap* _pred; + + typedef ExtendFindEnum TreeSet; + + IntNodeMap *_tree_set_index; + TreeSet *_tree_set; + + IntNodeMap *_delta1_index; + BinHeap *_delta1; + + IntNodeMap *_delta2_index; + BinHeap *_delta2; + + IntEdgeMap *_delta3_index; + BinHeap *_delta3; + + Value _delta_sum; + + void createStructures() { + _node_num = countNodes(_graph); + + if (!_matching) { + _matching = new MatchingMap(_graph); + } + if (!_node_potential) { + _node_potential = new NodePotential(_graph); + } + if (!_status) { + _status = new StatusMap(_graph); + } + if (!_pred) { + _pred = new PredMap(_graph); + } + if (!_tree_set) { + _tree_set_index = new IntNodeMap(_graph); + _tree_set = new TreeSet(*_tree_set_index); + } + if (!_delta1) { + _delta1_index = new IntNodeMap(_graph); + _delta1 = new BinHeap(*_delta1_index); + } + if (!_delta2) { + _delta2_index = new IntNodeMap(_graph); + _delta2 = new BinHeap(*_delta2_index); + } + if (!_delta3) { + _delta3_index = new IntEdgeMap(_graph); + _delta3 = new BinHeap(*_delta3_index); + } + } + + void destroyStructures() { + if (_matching) { + delete _matching; + } + if (_node_potential) { + delete _node_potential; + } + if (_status) { + delete _status; + } + if (_pred) { + delete _pred; + } + if (_tree_set) { + delete _tree_set_index; + delete _tree_set; + } + if (_delta1) { + delete _delta1_index; + delete _delta1; + } + if (_delta2) { + delete _delta2_index; + delete _delta2; + } + if (_delta3) { + delete _delta3_index; + delete _delta3; + } + } + + void matchedToEven(Node node, int tree) { + _tree_set->insert(node, tree); + _node_potential->set(node, (*_node_potential)[node] + _delta_sum); + _delta1->push(node, (*_node_potential)[node]); + + if (_delta2->state(node) == _delta2->IN_HEAP) { + _delta2->erase(node); + } + + for (InArcIt a(_graph, node); a != INVALID; ++a) { + Node v = _graph.source(a); + Value rw = (*_node_potential)[node] + (*_node_potential)[v] - + dualScale * _weight[a]; + if (node == v) { + if (_allow_loops && _graph.direction(a)) { + _delta3->push(a, rw / 2); + } + } else if ((*_status)[v] == EVEN) { + _delta3->push(a, rw / 2); + } else if ((*_status)[v] == MATCHED) { + if (_delta2->state(v) != _delta2->IN_HEAP) { + _pred->set(v, a); + _delta2->push(v, rw); + } else if ((*_delta2)[v] > rw) { + _pred->set(v, a); + _delta2->decrease(v, rw); + } + } + } + } + + void matchedToOdd(Node node, int tree) { + _tree_set->insert(node, tree); + _node_potential->set(node, (*_node_potential)[node] - _delta_sum); + + if (_delta2->state(node) == _delta2->IN_HEAP) { + _delta2->erase(node); + } + } + + void evenToMatched(Node node, int tree) { + _delta1->erase(node); + _node_potential->set(node, (*_node_potential)[node] - _delta_sum); + Arc min = INVALID; + Value minrw = std::numeric_limits::max(); + for (InArcIt a(_graph, node); a != INVALID; ++a) { + Node v = _graph.source(a); + Value rw = (*_node_potential)[node] + (*_node_potential)[v] - + dualScale * _weight[a]; + + if (node == v) { + if (_allow_loops && _graph.direction(a)) { + _delta3->erase(a); + } + } else if ((*_status)[v] == EVEN) { + _delta3->erase(a); + if (minrw > rw) { + min = _graph.oppositeArc(a); + minrw = rw; + } + } else if ((*_status)[v] == MATCHED) { + if ((*_pred)[v] == a) { + Arc mina = INVALID; + Value minrwa = std::numeric_limits::max(); + for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) { + Node va = _graph.target(aa); + if ((*_status)[va] != EVEN || + _tree_set->find(va) == tree) continue; + Value rwa = (*_node_potential)[v] + (*_node_potential)[va] - + dualScale * _weight[aa]; + if (minrwa > rwa) { + minrwa = rwa; + mina = aa; + } + } + if (mina != INVALID) { + _pred->set(v, mina); + _delta2->increase(v, minrwa); + } else { + _pred->set(v, INVALID); + _delta2->erase(v); + } + } + } + } + if (min != INVALID) { + _pred->set(node, min); + _delta2->push(node, minrw); + } else { + _pred->set(node, INVALID); + } + } + + void oddToMatched(Node node) { + _node_potential->set(node, (*_node_potential)[node] + _delta_sum); + Arc min = INVALID; + Value minrw = std::numeric_limits::max(); + for (InArcIt a(_graph, node); a != INVALID; ++a) { + Node v = _graph.source(a); + if ((*_status)[v] != EVEN) continue; + Value rw = (*_node_potential)[node] + (*_node_potential)[v] - + dualScale * _weight[a]; + + if (minrw > rw) { + min = _graph.oppositeArc(a); + minrw = rw; + } + } + if (min != INVALID) { + _pred->set(node, min); + _delta2->push(node, minrw); + } else { + _pred->set(node, INVALID); + } + } + + void alternatePath(Node even, int tree) { + Node odd; + + _status->set(even, MATCHED); + evenToMatched(even, tree); + + Arc prev = (*_matching)[even]; + while (prev != INVALID) { + odd = _graph.target(prev); + even = _graph.target((*_pred)[odd]); + _matching->set(odd, (*_pred)[odd]); + _status->set(odd, MATCHED); + oddToMatched(odd); + + prev = (*_matching)[even]; + _status->set(even, MATCHED); + _matching->set(even, _graph.oppositeArc((*_matching)[odd])); + evenToMatched(even, tree); + } + } + + void destroyTree(int tree) { + for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) { + if ((*_status)[n] == EVEN) { + _status->set(n, MATCHED); + evenToMatched(n, tree); + } else if ((*_status)[n] == ODD) { + _status->set(n, MATCHED); + oddToMatched(n); + } + } + _tree_set->eraseClass(tree); + } + + + void unmatchNode(const Node& node) { + int tree = _tree_set->find(node); + + alternatePath(node, tree); + destroyTree(tree); + + _matching->set(node, INVALID); + } + + + void augmentOnEdge(const Edge& edge) { + Node left = _graph.u(edge); + int left_tree = _tree_set->find(left); + + alternatePath(left, left_tree); + destroyTree(left_tree); + _matching->set(left, _graph.direct(edge, true)); + + Node right = _graph.v(edge); + int right_tree = _tree_set->find(right); + + alternatePath(right, right_tree); + destroyTree(right_tree); + _matching->set(right, _graph.direct(edge, false)); + } + + void augmentOnArc(const Arc& arc) { + Node left = _graph.source(arc); + _status->set(left, MATCHED); + _matching->set(left, arc); + _pred->set(left, arc); + + Node right = _graph.target(arc); + int right_tree = _tree_set->find(right); + + alternatePath(right, right_tree); + destroyTree(right_tree); + _matching->set(right, _graph.oppositeArc(arc)); + } + + void extendOnArc(const Arc& arc) { + Node base = _graph.target(arc); + int tree = _tree_set->find(base); + + Node odd = _graph.source(arc); + _tree_set->insert(odd, tree); + _status->set(odd, ODD); + matchedToOdd(odd, tree); + _pred->set(odd, arc); + + Node even = _graph.target((*_matching)[odd]); + _tree_set->insert(even, tree); + _status->set(even, EVEN); + matchedToEven(even, tree); + } + + void cycleOnEdge(const Edge& edge, int tree) { + Node nca = INVALID; + std::vector left_path, right_path; + + { + std::set left_set, right_set; + Node left = _graph.u(edge); + left_path.push_back(left); + left_set.insert(left); + + Node right = _graph.v(edge); + right_path.push_back(right); + right_set.insert(right); + + while (true) { + + if (left_set.find(right) != left_set.end()) { + nca = right; + break; + } + + if ((*_matching)[left] == INVALID) break; + + left = _graph.target((*_matching)[left]); + left_path.push_back(left); + left = _graph.target((*_pred)[left]); + left_path.push_back(left); + + left_set.insert(left); + + if (right_set.find(left) != right_set.end()) { + nca = left; + break; + } + + if ((*_matching)[right] == INVALID) break; + + right = _graph.target((*_matching)[right]); + right_path.push_back(right); + right = _graph.target((*_pred)[right]); + right_path.push_back(right); + + right_set.insert(right); + + } + + if (nca == INVALID) { + if ((*_matching)[left] == INVALID) { + nca = right; + while (left_set.find(nca) == left_set.end()) { + nca = _graph.target((*_matching)[nca]); + right_path.push_back(nca); + nca = _graph.target((*_pred)[nca]); + right_path.push_back(nca); + } + } else { + nca = left; + while (right_set.find(nca) == right_set.end()) { + nca = _graph.target((*_matching)[nca]); + left_path.push_back(nca); + nca = _graph.target((*_pred)[nca]); + left_path.push_back(nca); + } + } + } + } + + alternatePath(nca, tree); + Arc prev; + + prev = _graph.direct(edge, true); + for (int i = 0; left_path[i] != nca; i += 2) { + _matching->set(left_path[i], prev); + _status->set(left_path[i], MATCHED); + evenToMatched(left_path[i], tree); + + prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]); + _status->set(left_path[i + 1], MATCHED); + oddToMatched(left_path[i + 1]); + } + _matching->set(nca, prev); + + for (int i = 0; right_path[i] != nca; i += 2) { + _status->set(right_path[i], MATCHED); + evenToMatched(right_path[i], tree); + + _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]); + _status->set(right_path[i + 1], MATCHED); + oddToMatched(right_path[i + 1]); + } + + destroyTree(tree); + } + + void extractCycle(const Arc &arc) { + Node left = _graph.source(arc); + Node odd = _graph.target((*_matching)[left]); + Arc prev; + while (odd != left) { + Node even = _graph.target((*_matching)[odd]); + prev = (*_matching)[odd]; + odd = _graph.target((*_matching)[even]); + _matching->set(even, _graph.oppositeArc(prev)); + } + _matching->set(left, arc); + + Node right = _graph.target(arc); + int right_tree = _tree_set->find(right); + alternatePath(right, right_tree); + destroyTree(right_tree); + _matching->set(right, _graph.oppositeArc(arc)); + } + + public: + + /// \brief Constructor + /// + /// Constructor. + MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight, + bool allow_loops = true) + : _graph(graph), _weight(weight), _matching(0), + _node_potential(0), _node_num(0), _allow_loops(allow_loops), + _status(0), _pred(0), + _tree_set_index(0), _tree_set(0), + + _delta1_index(0), _delta1(0), + _delta2_index(0), _delta2(0), + _delta3_index(0), _delta3(0), + + _delta_sum() {} + + ~MaxWeightedFractionalMatching() { + destroyStructures(); + } + + /// \name Execution Control + /// The simplest way to execute the algorithm is to use the + /// \ref run() member function. + + ///@{ + + /// \brief Initialize the algorithm + /// + /// This function initializes the algorithm. + void init() { + createStructures(); + + for (NodeIt n(_graph); n != INVALID; ++n) { + (*_delta1_index)[n] = _delta1->PRE_HEAP; + (*_delta2_index)[n] = _delta2->PRE_HEAP; + } + for (EdgeIt e(_graph); e != INVALID; ++e) { + (*_delta3_index)[e] = _delta3->PRE_HEAP; + } + + for (NodeIt n(_graph); n != INVALID; ++n) { + Value max = 0; + for (OutArcIt e(_graph, n); e != INVALID; ++e) { + if (_graph.target(e) == n && !_allow_loops) continue; + if ((dualScale * _weight[e]) / 2 > max) { + max = (dualScale * _weight[e]) / 2; + } + } + _node_potential->set(n, max); + _delta1->push(n, max); + + _tree_set->insert(n); + + _matching->set(n, INVALID); + _status->set(n, EVEN); + } + + for (EdgeIt e(_graph); e != INVALID; ++e) { + Node left = _graph.u(e); + Node right = _graph.v(e); + if (left == right && !_allow_loops) continue; + _delta3->push(e, ((*_node_potential)[left] + + (*_node_potential)[right] - + dualScale * _weight[e]) / 2); + } + } + + /// \brief Start the algorithm + /// + /// This function starts the algorithm. + /// + /// \pre \ref init() must be called before using this function. + void start() { + enum OpType { + D1, D2, D3 + }; + + int unmatched = _node_num; + while (unmatched > 0) { + Value d1 = !_delta1->empty() ? + _delta1->prio() : std::numeric_limits::max(); + + Value d2 = !_delta2->empty() ? + _delta2->prio() : std::numeric_limits::max(); + + Value d3 = !_delta3->empty() ? + _delta3->prio() : std::numeric_limits::max(); + + _delta_sum = d3; OpType ot = D3; + if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; } + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } + + switch (ot) { + case D1: + { + Node n = _delta1->top(); + unmatchNode(n); + --unmatched; + } + break; + case D2: + { + Node n = _delta2->top(); + Arc a = (*_pred)[n]; + if ((*_matching)[n] == INVALID) { + augmentOnArc(a); + --unmatched; + } else { + Node v = _graph.target((*_matching)[n]); + if ((*_matching)[n] != + _graph.oppositeArc((*_matching)[v])) { + extractCycle(a); + --unmatched; + } else { + extendOnArc(a); + } + } + } break; + case D3: + { + Edge e = _delta3->top(); + + Node left = _graph.u(e); + Node right = _graph.v(e); + + int left_tree = _tree_set->find(left); + int right_tree = _tree_set->find(right); + + if (left_tree == right_tree) { + cycleOnEdge(e, left_tree); + --unmatched; + } else { + augmentOnEdge(e); + unmatched -= 2; + } + } break; + } + } + } + + /// \brief Run the algorithm. + /// + /// This method runs the \c %MaxWeightedMatching algorithm. + /// + /// \note mwfm.run() is just a shortcut of the following code. + /// \code + /// mwfm.init(); + /// mwfm.start(); + /// \endcode + void run() { + init(); + start(); + } + + /// @} + + /// \name Primal Solution + /// Functions to get the primal solution, i.e. the maximum weighted + /// matching.\n + /// Either \ref run() or \ref start() function should be called before + /// using them. + + /// @{ + + /// \brief Return the weight of the matching. + /// + /// This function returns the weight of the found matching. This + /// value is scaled by \ref primalScale "primal scale". + /// + /// \pre Either run() or start() must be called before using this function. + Value matchingWeight() const { + Value sum = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_matching)[n] != INVALID) { + sum += _weight[(*_matching)[n]]; + } + } + return sum * primalScale / 2; + } + + /// \brief Return the number of covered nodes in the matching. + /// + /// This function returns the number of covered nodes in the matching. + /// + /// \pre Either run() or start() must be called before using this function. + int matchingSize() const { + int num = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_matching)[n] != INVALID) { + ++num; + } + } + return num; + } + + /// \brief Return \c true if the given edge is in the matching. + /// + /// This function returns \c true if the given edge is in the + /// found matching. The result is scaled by \ref primalScale + /// "primal scale". + /// + /// \pre Either run() or start() must be called before using this function. + Value matching(const Edge& edge) const { + return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0) + * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0) + * primalScale / 2; + } + + /// \brief Return the fractional matching arc (or edge) incident + /// to the given node. + /// + /// This function returns one of the fractional matching arc (or + /// edge) incident to the given node in the found matching or \c + /// INVALID if the node is not covered by the matching or if the + /// node is on an odd length cycle then it is the successor edge + /// on the cycle. + /// + /// \pre Either run() or start() must be called before using this function. + Arc matching(const Node& node) const { + return (*_matching)[node]; + } + + /// \brief Return a const reference to the matching map. + /// + /// This function returns a const reference to a node map that stores + /// the matching arc (or edge) incident to each node. + const MatchingMap& matchingMap() const { + return *_matching; + } + + /// @} + + /// \name Dual Solution + /// Functions to get the dual solution.\n + /// Either \ref run() or \ref start() function should be called before + /// using them. + + /// @{ + + /// \brief Return the value of the dual solution. + /// + /// This function returns the value of the dual solution. + /// It should be equal to the primal value scaled by \ref dualScale + /// "dual scale". + /// + /// \pre Either run() or start() must be called before using this function. + Value dualValue() const { + Value sum = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + sum += nodeValue(n); + } + return sum; + } + + /// \brief Return the dual value (potential) of the given node. + /// + /// This function returns the dual value (potential) of the given node. + /// + /// \pre Either run() or start() must be called before using this function. + Value nodeValue(const Node& n) const { + return (*_node_potential)[n]; + } + + /// @} + + }; + + /// \ingroup matching + /// + /// \brief Weighted fractional perfect matching in general graphs + /// + /// This class provides an efficient implementation of fractional + /// matching algorithm. The implementation is based on extensive use + /// of priority queues and provides \f$O(nm\log n)\f$ time + /// complexity. + /// + /// The maximum weighted fractional perfect matching is a relaxation + /// of the maximum weighted perfect matching problem where the odd + /// set constraints are omitted. + /// It can be formulated with the following linear program. + /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f] + /// \f[x_e \ge 0\quad \forall e\in E\f] + /// \f[\max \sum_{e\in E}x_ew_e\f] + /// where \f$\delta(X)\f$ is the set of edges incident to a node in + /// \f$X\f$. The result must be the union of a matching with one + /// value edges and a set of odd length cycles with half value edges. + /// + /// The algorithm calculates an optimal fractional matching and a + /// proof of the optimality. The solution of the dual problem can be + /// used to check the result of the algorithm. The dual linear + /// problem is the following. + /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f] + /// \f[\min \sum_{u \in V}y_u \f] /// + /// + /// The algorithm can be executed with the run() function. + /// After it the matching (the primal solution) and the dual solution + /// can be obtained using the query functions. + + /// If the value type is integer, then the primal and the dual + /// solutions are multiplied by + /// \ref MaxWeightedMatching::primalScale "2" and + /// \ref MaxWeightedMatching::dualScale "4" respectively. + /// + /// \tparam GR The undirected graph type the algorithm runs on. + /// \tparam WM The type edge weight map. The default type is + /// \ref concepts::Graph::EdgeMap "GR::EdgeMap". +#ifdef DOXYGEN + template +#else + template > +#endif + class MaxWeightedPerfectFractionalMatching { + public: + + /// The graph type of the algorithm + typedef GR Graph; + /// The type of the edge weight map + typedef WM WeightMap; + /// The value type of the edge weights + typedef typename WeightMap::Value Value; + + /// The type of the matching map + typedef typename Graph::template NodeMap + MatchingMap; + + /// \brief Scaling factor for primal solution + /// + /// Scaling factor for primal solution. It is equal to 2 or 1 + /// according to the value type. + static const int primalScale = + std::numeric_limits::is_integer ? 2 : 1; + + /// \brief Scaling factor for dual solution + /// + /// Scaling factor for dual solution. It is equal to 4 or 1 + /// according to the value type. + static const int dualScale = + std::numeric_limits::is_integer ? 4 : 1; + + private: + + TEMPLATE_GRAPH_TYPEDEFS(Graph); + + typedef typename Graph::template NodeMap NodePotential; + + const Graph& _graph; + const WeightMap& _weight; + + MatchingMap* _matching; + NodePotential* _node_potential; + + int _node_num; + bool _allow_loops; + + enum Status { + EVEN = -1, MATCHED = 0, ODD = 1 + }; + + typedef typename Graph::template NodeMap StatusMap; + StatusMap* _status; + + typedef typename Graph::template NodeMap PredMap; + PredMap* _pred; + + typedef ExtendFindEnum TreeSet; + + IntNodeMap *_tree_set_index; + TreeSet *_tree_set; + + IntNodeMap *_delta2_index; + BinHeap *_delta2; + + IntEdgeMap *_delta3_index; + BinHeap *_delta3; + + Value _delta_sum; + + void createStructures() { + _node_num = countNodes(_graph); + + if (!_matching) { + _matching = new MatchingMap(_graph); + } + if (!_node_potential) { + _node_potential = new NodePotential(_graph); + } + if (!_status) { + _status = new StatusMap(_graph); + } + if (!_pred) { + _pred = new PredMap(_graph); + } + if (!_tree_set) { + _tree_set_index = new IntNodeMap(_graph); + _tree_set = new TreeSet(*_tree_set_index); + } + if (!_delta2) { + _delta2_index = new IntNodeMap(_graph); + _delta2 = new BinHeap(*_delta2_index); + } + if (!_delta3) { + _delta3_index = new IntEdgeMap(_graph); + _delta3 = new BinHeap(*_delta3_index); + } + } + + void destroyStructures() { + if (_matching) { + delete _matching; + } + if (_node_potential) { + delete _node_potential; + } + if (_status) { + delete _status; + } + if (_pred) { + delete _pred; + } + if (_tree_set) { + delete _tree_set_index; + delete _tree_set; + } + if (_delta2) { + delete _delta2_index; + delete _delta2; + } + if (_delta3) { + delete _delta3_index; + delete _delta3; + } + } + + void matchedToEven(Node node, int tree) { + _tree_set->insert(node, tree); + _node_potential->set(node, (*_node_potential)[node] + _delta_sum); + + if (_delta2->state(node) == _delta2->IN_HEAP) { + _delta2->erase(node); + } + + for (InArcIt a(_graph, node); a != INVALID; ++a) { + Node v = _graph.source(a); + Value rw = (*_node_potential)[node] + (*_node_potential)[v] - + dualScale * _weight[a]; + if (node == v) { + if (_allow_loops && _graph.direction(a)) { + _delta3->push(a, rw / 2); + } + } else if ((*_status)[v] == EVEN) { + _delta3->push(a, rw / 2); + } else if ((*_status)[v] == MATCHED) { + if (_delta2->state(v) != _delta2->IN_HEAP) { + _pred->set(v, a); + _delta2->push(v, rw); + } else if ((*_delta2)[v] > rw) { + _pred->set(v, a); + _delta2->decrease(v, rw); + } + } + } + } + + void matchedToOdd(Node node, int tree) { + _tree_set->insert(node, tree); + _node_potential->set(node, (*_node_potential)[node] - _delta_sum); + + if (_delta2->state(node) == _delta2->IN_HEAP) { + _delta2->erase(node); + } + } + + void evenToMatched(Node node, int tree) { + _node_potential->set(node, (*_node_potential)[node] - _delta_sum); + Arc min = INVALID; + Value minrw = std::numeric_limits::max(); + for (InArcIt a(_graph, node); a != INVALID; ++a) { + Node v = _graph.source(a); + Value rw = (*_node_potential)[node] + (*_node_potential)[v] - + dualScale * _weight[a]; + + if (node == v) { + if (_allow_loops && _graph.direction(a)) { + _delta3->erase(a); + } + } else if ((*_status)[v] == EVEN) { + _delta3->erase(a); + if (minrw > rw) { + min = _graph.oppositeArc(a); + minrw = rw; + } + } else if ((*_status)[v] == MATCHED) { + if ((*_pred)[v] == a) { + Arc mina = INVALID; + Value minrwa = std::numeric_limits::max(); + for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) { + Node va = _graph.target(aa); + if ((*_status)[va] != EVEN || + _tree_set->find(va) == tree) continue; + Value rwa = (*_node_potential)[v] + (*_node_potential)[va] - + dualScale * _weight[aa]; + if (minrwa > rwa) { + minrwa = rwa; + mina = aa; + } + } + if (mina != INVALID) { + _pred->set(v, mina); + _delta2->increase(v, minrwa); + } else { + _pred->set(v, INVALID); + _delta2->erase(v); + } + } + } + } + if (min != INVALID) { + _pred->set(node, min); + _delta2->push(node, minrw); + } else { + _pred->set(node, INVALID); + } + } + + void oddToMatched(Node node) { + _node_potential->set(node, (*_node_potential)[node] + _delta_sum); + Arc min = INVALID; + Value minrw = std::numeric_limits::max(); + for (InArcIt a(_graph, node); a != INVALID; ++a) { + Node v = _graph.source(a); + if ((*_status)[v] != EVEN) continue; + Value rw = (*_node_potential)[node] + (*_node_potential)[v] - + dualScale * _weight[a]; + + if (minrw > rw) { + min = _graph.oppositeArc(a); + minrw = rw; + } + } + if (min != INVALID) { + _pred->set(node, min); + _delta2->push(node, minrw); + } else { + _pred->set(node, INVALID); + } + } + + void alternatePath(Node even, int tree) { + Node odd; + + _status->set(even, MATCHED); + evenToMatched(even, tree); + + Arc prev = (*_matching)[even]; + while (prev != INVALID) { + odd = _graph.target(prev); + even = _graph.target((*_pred)[odd]); + _matching->set(odd, (*_pred)[odd]); + _status->set(odd, MATCHED); + oddToMatched(odd); + + prev = (*_matching)[even]; + _status->set(even, MATCHED); + _matching->set(even, _graph.oppositeArc((*_matching)[odd])); + evenToMatched(even, tree); + } + } + + void destroyTree(int tree) { + for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) { + if ((*_status)[n] == EVEN) { + _status->set(n, MATCHED); + evenToMatched(n, tree); + } else if ((*_status)[n] == ODD) { + _status->set(n, MATCHED); + oddToMatched(n); + } + } + _tree_set->eraseClass(tree); + } + + void augmentOnEdge(const Edge& edge) { + Node left = _graph.u(edge); + int left_tree = _tree_set->find(left); + + alternatePath(left, left_tree); + destroyTree(left_tree); + _matching->set(left, _graph.direct(edge, true)); + + Node right = _graph.v(edge); + int right_tree = _tree_set->find(right); + + alternatePath(right, right_tree); + destroyTree(right_tree); + _matching->set(right, _graph.direct(edge, false)); + } + + void augmentOnArc(const Arc& arc) { + Node left = _graph.source(arc); + _status->set(left, MATCHED); + _matching->set(left, arc); + _pred->set(left, arc); + + Node right = _graph.target(arc); + int right_tree = _tree_set->find(right); + + alternatePath(right, right_tree); + destroyTree(right_tree); + _matching->set(right, _graph.oppositeArc(arc)); + } + + void extendOnArc(const Arc& arc) { + Node base = _graph.target(arc); + int tree = _tree_set->find(base); + + Node odd = _graph.source(arc); + _tree_set->insert(odd, tree); + _status->set(odd, ODD); + matchedToOdd(odd, tree); + _pred->set(odd, arc); + + Node even = _graph.target((*_matching)[odd]); + _tree_set->insert(even, tree); + _status->set(even, EVEN); + matchedToEven(even, tree); + } + + void cycleOnEdge(const Edge& edge, int tree) { + Node nca = INVALID; + std::vector left_path, right_path; + + { + std::set left_set, right_set; + Node left = _graph.u(edge); + left_path.push_back(left); + left_set.insert(left); + + Node right = _graph.v(edge); + right_path.push_back(right); + right_set.insert(right); + + while (true) { + + if (left_set.find(right) != left_set.end()) { + nca = right; + break; + } + + if ((*_matching)[left] == INVALID) break; + + left = _graph.target((*_matching)[left]); + left_path.push_back(left); + left = _graph.target((*_pred)[left]); + left_path.push_back(left); + + left_set.insert(left); + + if (right_set.find(left) != right_set.end()) { + nca = left; + break; + } + + if ((*_matching)[right] == INVALID) break; + + right = _graph.target((*_matching)[right]); + right_path.push_back(right); + right = _graph.target((*_pred)[right]); + right_path.push_back(right); + + right_set.insert(right); + + } + + if (nca == INVALID) { + if ((*_matching)[left] == INVALID) { + nca = right; + while (left_set.find(nca) == left_set.end()) { + nca = _graph.target((*_matching)[nca]); + right_path.push_back(nca); + nca = _graph.target((*_pred)[nca]); + right_path.push_back(nca); + } + } else { + nca = left; + while (right_set.find(nca) == right_set.end()) { + nca = _graph.target((*_matching)[nca]); + left_path.push_back(nca); + nca = _graph.target((*_pred)[nca]); + left_path.push_back(nca); + } + } + } + } + + alternatePath(nca, tree); + Arc prev; + + prev = _graph.direct(edge, true); + for (int i = 0; left_path[i] != nca; i += 2) { + _matching->set(left_path[i], prev); + _status->set(left_path[i], MATCHED); + evenToMatched(left_path[i], tree); + + prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]); + _status->set(left_path[i + 1], MATCHED); + oddToMatched(left_path[i + 1]); + } + _matching->set(nca, prev); + + for (int i = 0; right_path[i] != nca; i += 2) { + _status->set(right_path[i], MATCHED); + evenToMatched(right_path[i], tree); + + _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]); + _status->set(right_path[i + 1], MATCHED); + oddToMatched(right_path[i + 1]); + } + + destroyTree(tree); + } + + void extractCycle(const Arc &arc) { + Node left = _graph.source(arc); + Node odd = _graph.target((*_matching)[left]); + Arc prev; + while (odd != left) { + Node even = _graph.target((*_matching)[odd]); + prev = (*_matching)[odd]; + odd = _graph.target((*_matching)[even]); + _matching->set(even, _graph.oppositeArc(prev)); + } + _matching->set(left, arc); + + Node right = _graph.target(arc); + int right_tree = _tree_set->find(right); + alternatePath(right, right_tree); + destroyTree(right_tree); + _matching->set(right, _graph.oppositeArc(arc)); + } + + public: + + /// \brief Constructor + /// + /// Constructor. + MaxWeightedPerfectFractionalMatching(const Graph& graph, + const WeightMap& weight, + bool allow_loops = true) + : _graph(graph), _weight(weight), _matching(0), + _node_potential(0), _node_num(0), _allow_loops(allow_loops), + _status(0), _pred(0), + _tree_set_index(0), _tree_set(0), + + _delta2_index(0), _delta2(0), + _delta3_index(0), _delta3(0), + + _delta_sum() {} + + ~MaxWeightedPerfectFractionalMatching() { + destroyStructures(); + } + + /// \name Execution Control + /// The simplest way to execute the algorithm is to use the + /// \ref run() member function. + + ///@{ + + /// \brief Initialize the algorithm + /// + /// This function initializes the algorithm. + void init() { + createStructures(); + + for (NodeIt n(_graph); n != INVALID; ++n) { + (*_delta2_index)[n] = _delta2->PRE_HEAP; + } + for (EdgeIt e(_graph); e != INVALID; ++e) { + (*_delta3_index)[e] = _delta3->PRE_HEAP; + } + + for (NodeIt n(_graph); n != INVALID; ++n) { + Value max = - std::numeric_limits::max(); + for (OutArcIt e(_graph, n); e != INVALID; ++e) { + if (_graph.target(e) == n && !_allow_loops) continue; + if ((dualScale * _weight[e]) / 2 > max) { + max = (dualScale * _weight[e]) / 2; + } + } + _node_potential->set(n, max); + + _tree_set->insert(n); + + _matching->set(n, INVALID); + _status->set(n, EVEN); + } + + for (EdgeIt e(_graph); e != INVALID; ++e) { + Node left = _graph.u(e); + Node right = _graph.v(e); + if (left == right && !_allow_loops) continue; + _delta3->push(e, ((*_node_potential)[left] + + (*_node_potential)[right] - + dualScale * _weight[e]) / 2); + } + } + + /// \brief Start the algorithm + /// + /// This function starts the algorithm. + /// + /// \pre \ref init() must be called before using this function. + bool start() { + enum OpType { + D2, D3 + }; + + int unmatched = _node_num; + while (unmatched > 0) { + Value d2 = !_delta2->empty() ? + _delta2->prio() : std::numeric_limits::max(); + + Value d3 = !_delta3->empty() ? + _delta3->prio() : std::numeric_limits::max(); + + _delta_sum = d3; OpType ot = D3; + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } + + if (_delta_sum == std::numeric_limits::max()) { + return false; + } + + switch (ot) { + case D2: + { + Node n = _delta2->top(); + Arc a = (*_pred)[n]; + if ((*_matching)[n] == INVALID) { + augmentOnArc(a); + --unmatched; + } else { + Node v = _graph.target((*_matching)[n]); + if ((*_matching)[n] != + _graph.oppositeArc((*_matching)[v])) { + extractCycle(a); + --unmatched; + } else { + extendOnArc(a); + } + } + } break; + case D3: + { + Edge e = _delta3->top(); + + Node left = _graph.u(e); + Node right = _graph.v(e); + + int left_tree = _tree_set->find(left); + int right_tree = _tree_set->find(right); + + if (left_tree == right_tree) { + cycleOnEdge(e, left_tree); + --unmatched; + } else { + augmentOnEdge(e); + unmatched -= 2; + } + } break; + } + } + return true; + } + + /// \brief Run the algorithm. + /// + /// This method runs the \c %MaxWeightedMatching algorithm. + /// + /// \note mwfm.run() is just a shortcut of the following code. + /// \code + /// mwpfm.init(); + /// mwpfm.start(); + /// \endcode + bool run() { + init(); + return start(); + } + + /// @} + + /// \name Primal Solution + /// Functions to get the primal solution, i.e. the maximum weighted + /// matching.\n + /// Either \ref run() or \ref start() function should be called before + /// using them. + + /// @{ + + /// \brief Return the weight of the matching. + /// + /// This function returns the weight of the found matching. This + /// value is scaled by \ref primalScale "primal scale". + /// + /// \pre Either run() or start() must be called before using this function. + Value matchingWeight() const { + Value sum = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_matching)[n] != INVALID) { + sum += _weight[(*_matching)[n]]; + } + } + return sum * primalScale / 2; + } + + /// \brief Return the number of covered nodes in the matching. + /// + /// This function returns the number of covered nodes in the matching. + /// + /// \pre Either run() or start() must be called before using this function. + int matchingSize() const { + int num = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + if ((*_matching)[n] != INVALID) { + ++num; + } + } + return num; + } + + /// \brief Return \c true if the given edge is in the matching. + /// + /// This function returns \c true if the given edge is in the + /// found matching. The result is scaled by \ref primalScale + /// "primal scale". + /// + /// \pre Either run() or start() must be called before using this function. + Value matching(const Edge& edge) const { + return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0) + * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0) + * primalScale / 2; + } + + /// \brief Return the fractional matching arc (or edge) incident + /// to the given node. + /// + /// This function returns one of the fractional matching arc (or + /// edge) incident to the given node in the found matching or \c + /// INVALID if the node is not covered by the matching or if the + /// node is on an odd length cycle then it is the successor edge + /// on the cycle. + /// + /// \pre Either run() or start() must be called before using this function. + Arc matching(const Node& node) const { + return (*_matching)[node]; + } + + /// \brief Return a const reference to the matching map. + /// + /// This function returns a const reference to a node map that stores + /// the matching arc (or edge) incident to each node. + const MatchingMap& matchingMap() const { + return *_matching; + } + + /// @} + + /// \name Dual Solution + /// Functions to get the dual solution.\n + /// Either \ref run() or \ref start() function should be called before + /// using them. + + /// @{ + + /// \brief Return the value of the dual solution. + /// + /// This function returns the value of the dual solution. + /// It should be equal to the primal value scaled by \ref dualScale + /// "dual scale". + /// + /// \pre Either run() or start() must be called before using this function. + Value dualValue() const { + Value sum = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + sum += nodeValue(n); + } + return sum; + } + + /// \brief Return the dual value (potential) of the given node. + /// + /// This function returns the dual value (potential) of the given node. + /// + /// \pre Either run() or start() must be called before using this function. + Value nodeValue(const Node& n) const { + return (*_node_potential)[n]; + } + + /// @} + + }; + +} //END OF NAMESPACE LEMON + +#endif //LEMON_FRACTIONAL_MATCHING_H diff -r 0513ccfea967 -r 636dadefe1e6 test/CMakeLists.txt --- a/test/CMakeLists.txt Sun Sep 20 21:38:24 2009 +0200 +++ b/test/CMakeLists.txt Fri Sep 25 21:51:36 2009 +0200 @@ -21,6 +21,7 @@ edge_set_test error_test euler_test + fractional_matching_test gomory_hu_test graph_copy_test graph_test diff -r 0513ccfea967 -r 636dadefe1e6 test/Makefile.am --- a/test/Makefile.am Sun Sep 20 21:38:24 2009 +0200 +++ b/test/Makefile.am Fri Sep 25 21:51:36 2009 +0200 @@ -19,6 +19,7 @@ test/edge_set_test \ test/error_test \ test/euler_test \ + test/fractional_matching_test \ test/gomory_hu_test \ test/graph_copy_test \ test/graph_test \ @@ -65,6 +66,7 @@ test_edge_set_test_SOURCES = test/edge_set_test.cc test_error_test_SOURCES = test/error_test.cc test_euler_test_SOURCES = test/euler_test.cc +test_fractional_matching_test_SOURCES = test/fractional_matching_test.cc test_gomory_hu_test_SOURCES = test/gomory_hu_test.cc test_graph_copy_test_SOURCES = test/graph_copy_test.cc test_graph_test_SOURCES = test/graph_test.cc diff -r 0513ccfea967 -r 636dadefe1e6 test/fractional_matching_test.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/fractional_matching_test.cc Fri Sep 25 21:51:36 2009 +0200 @@ -0,0 +1,502 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2009 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "test_tools.h" + +using namespace std; +using namespace lemon; + +GRAPH_TYPEDEFS(SmartGraph); + + +const int lgfn = 4; +const std::string lgf[lgfn] = { + "@nodes\n" + "label\n" + "0\n" + "1\n" + "2\n" + "3\n" + "4\n" + "5\n" + "6\n" + "7\n" + "@edges\n" + " label weight\n" + "7 4 0 984\n" + "0 7 1 73\n" + "7 1 2 204\n" + "2 3 3 583\n" + "2 7 4 565\n" + "2 1 5 582\n" + "0 4 6 551\n" + "2 5 7 385\n" + "1 5 8 561\n" + "5 3 9 484\n" + "7 5 10 904\n" + "3 6 11 47\n" + "7 6 12 888\n" + "3 0 13 747\n" + "6 1 14 310\n", + + "@nodes\n" + "label\n" + "0\n" + "1\n" + "2\n" + "3\n" + "4\n" + "5\n" + "6\n" + "7\n" + "@edges\n" + " label weight\n" + "2 5 0 710\n" + "0 5 1 241\n" + "2 4 2 856\n" + "2 6 3 762\n" + "4 1 4 747\n" + "6 1 5 962\n" + "4 7 6 723\n" + "1 7 7 661\n" + "2 3 8 376\n" + "1 0 9 416\n" + "6 7 10 391\n", + + "@nodes\n" + "label\n" + "0\n" + "1\n" + "2\n" + "3\n" + "4\n" + "5\n" + "6\n" + "7\n" + "@edges\n" + " label weight\n" + "6 2 0 553\n" + "0 7 1 653\n" + "6 3 2 22\n" + "4 7 3 846\n" + "7 2 4 981\n" + "7 6 5 250\n" + "5 2 6 539\n", + + "@nodes\n" + "label\n" + "0\n" + "@edges\n" + " label weight\n" + "0 0 0 100\n" +}; + +void checkMaxFractionalMatchingCompile() +{ + typedef concepts::Graph Graph; + typedef Graph::Node Node; + typedef Graph::Edge Edge; + + Graph g; + Node n; + Edge e; + + MaxFractionalMatching mat_test(g); + const MaxFractionalMatching& + const_mat_test = mat_test; + + mat_test.init(); + mat_test.start(); + mat_test.start(true); + mat_test.startPerfect(); + mat_test.startPerfect(true); + mat_test.run(); + mat_test.run(true); + mat_test.runPerfect(); + mat_test.runPerfect(true); + + const_mat_test.matchingSize(); + const_mat_test.matching(e); + const_mat_test.matching(n); + const MaxFractionalMatching::MatchingMap& mmap = + const_mat_test.matchingMap(); + e = mmap[n]; + + const_mat_test.barrier(n); +} + +void checkMaxWeightedFractionalMatchingCompile() +{ + typedef concepts::Graph Graph; + typedef Graph::Node Node; + typedef Graph::Edge Edge; + typedef Graph::EdgeMap WeightMap; + + Graph g; + Node n; + Edge e; + WeightMap w(g); + + MaxWeightedFractionalMatching mat_test(g, w); + const MaxWeightedFractionalMatching& + const_mat_test = mat_test; + + mat_test.init(); + mat_test.start(); + mat_test.run(); + + const_mat_test.matchingWeight(); + const_mat_test.matchingSize(); + const_mat_test.matching(e); + const_mat_test.matching(n); + const MaxWeightedFractionalMatching::MatchingMap& mmap = + const_mat_test.matchingMap(); + e = mmap[n]; + + const_mat_test.dualValue(); + const_mat_test.nodeValue(n); +} + +void checkMaxWeightedPerfectFractionalMatchingCompile() +{ + typedef concepts::Graph Graph; + typedef Graph::Node Node; + typedef Graph::Edge Edge; + typedef Graph::EdgeMap WeightMap; + + Graph g; + Node n; + Edge e; + WeightMap w(g); + + MaxWeightedPerfectFractionalMatching mat_test(g, w); + const MaxWeightedPerfectFractionalMatching& + const_mat_test = mat_test; + + mat_test.init(); + mat_test.start(); + mat_test.run(); + + const_mat_test.matchingWeight(); + const_mat_test.matching(e); + const_mat_test.matching(n); + const MaxWeightedPerfectFractionalMatching::MatchingMap& mmap = + const_mat_test.matchingMap(); + e = mmap[n]; + + const_mat_test.dualValue(); + const_mat_test.nodeValue(n); +} + +void checkFractionalMatching(const SmartGraph& graph, + const MaxFractionalMatching& mfm, + bool allow_loops = true) { + int pv = 0; + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + int indeg = 0; + for (InArcIt a(graph, n); a != INVALID; ++a) { + if (mfm.matching(graph.source(a)) == a) { + ++indeg; + } + } + if (mfm.matching(n) != INVALID) { + check(indeg == 1, "Invalid matching"); + ++pv; + } else { + check(indeg == 0, "Invalid matching"); + } + } + check(pv == mfm.matchingSize(), "Wrong matching size"); + + SmartGraph::NodeMap processed(graph, false); + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + if (processed[n]) continue; + processed[n] = true; + if (mfm.matching(n) == INVALID) continue; + int num = 1; + Node v = graph.target(mfm.matching(n)); + while (v != n) { + processed[v] = true; + ++num; + v = graph.target(mfm.matching(v)); + } + check(num == 2 || num % 2 == 1, "Wrong cycle size"); + check(allow_loops || num != 1, "Wrong cycle size"); + } + + int anum = 0, bnum = 0; + SmartGraph::NodeMap neighbours(graph, false); + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + if (!mfm.barrier(n)) continue; + ++anum; + for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) { + Node u = graph.source(a); + if (!allow_loops && u == n) continue; + if (!neighbours[u]) { + neighbours[u] = true; + ++bnum; + } + } + } + check(anum - bnum + mfm.matchingSize() == countNodes(graph), + "Wrong barrier"); +} + +void checkPerfectFractionalMatching(const SmartGraph& graph, + const MaxFractionalMatching& mfm, + bool perfect, bool allow_loops = true) { + if (perfect) { + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + int indeg = 0; + for (InArcIt a(graph, n); a != INVALID; ++a) { + if (mfm.matching(graph.source(a)) == a) { + ++indeg; + } + } + check(mfm.matching(n) != INVALID, "Invalid matching"); + check(indeg == 1, "Invalid matching"); + } + } else { + int anum = 0, bnum = 0; + SmartGraph::NodeMap neighbours(graph, false); + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + if (!mfm.barrier(n)) continue; + ++anum; + for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) { + Node u = graph.source(a); + if (!allow_loops && u == n) continue; + if (!neighbours[u]) { + neighbours[u] = true; + ++bnum; + } + } + } + check(anum - bnum > 0, "Wrong barrier"); + } +} + +void checkWeightedFractionalMatching(const SmartGraph& graph, + const SmartGraph::EdgeMap& weight, + const MaxWeightedFractionalMatching& mwfm, + bool allow_loops = true) { + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) { + if (graph.u(e) == graph.v(e) && !allow_loops) continue; + int rw = mwfm.nodeValue(graph.u(e)) + mwfm.nodeValue(graph.v(e)) + - weight[e] * mwfm.dualScale; + + check(rw >= 0, "Negative reduced weight"); + check(rw == 0 || !mwfm.matching(e), + "Non-zero reduced weight on matching edge"); + } + + int pv = 0; + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + int indeg = 0; + for (InArcIt a(graph, n); a != INVALID; ++a) { + if (mwfm.matching(graph.source(a)) == a) { + ++indeg; + } + } + check(indeg <= 1, "Invalid matching"); + if (mwfm.matching(n) != INVALID) { + check(mwfm.nodeValue(n) >= 0, "Invalid node value"); + check(indeg == 1, "Invalid matching"); + pv += weight[mwfm.matching(n)]; + SmartGraph::Node o = graph.target(mwfm.matching(n)); + } else { + check(mwfm.nodeValue(n) == 0, "Invalid matching"); + check(indeg == 0, "Invalid matching"); + } + } + + int dv = 0; + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + dv += mwfm.nodeValue(n); + } + + check(pv * mwfm.dualScale == dv * 2, "Wrong duality"); + + SmartGraph::NodeMap processed(graph, false); + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + if (processed[n]) continue; + processed[n] = true; + if (mwfm.matching(n) == INVALID) continue; + int num = 1; + Node v = graph.target(mwfm.matching(n)); + while (v != n) { + processed[v] = true; + ++num; + v = graph.target(mwfm.matching(v)); + } + check(num == 2 || num % 2 == 1, "Wrong cycle size"); + check(allow_loops || num != 1, "Wrong cycle size"); + } + + return; +} + +void checkWeightedPerfectFractionalMatching(const SmartGraph& graph, + const SmartGraph::EdgeMap& weight, + const MaxWeightedPerfectFractionalMatching& mwpfm, + bool allow_loops = true) { + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) { + if (graph.u(e) == graph.v(e) && !allow_loops) continue; + int rw = mwpfm.nodeValue(graph.u(e)) + mwpfm.nodeValue(graph.v(e)) + - weight[e] * mwpfm.dualScale; + + check(rw >= 0, "Negative reduced weight"); + check(rw == 0 || !mwpfm.matching(e), + "Non-zero reduced weight on matching edge"); + } + + int pv = 0; + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + int indeg = 0; + for (InArcIt a(graph, n); a != INVALID; ++a) { + if (mwpfm.matching(graph.source(a)) == a) { + ++indeg; + } + } + check(mwpfm.matching(n) != INVALID, "Invalid perfect matching"); + check(indeg == 1, "Invalid perfect matching"); + pv += weight[mwpfm.matching(n)]; + SmartGraph::Node o = graph.target(mwpfm.matching(n)); + } + + int dv = 0; + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + dv += mwpfm.nodeValue(n); + } + + check(pv * mwpfm.dualScale == dv * 2, "Wrong duality"); + + SmartGraph::NodeMap processed(graph, false); + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) { + if (processed[n]) continue; + processed[n] = true; + if (mwpfm.matching(n) == INVALID) continue; + int num = 1; + Node v = graph.target(mwpfm.matching(n)); + while (v != n) { + processed[v] = true; + ++num; + v = graph.target(mwpfm.matching(v)); + } + check(num == 2 || num % 2 == 1, "Wrong cycle size"); + check(allow_loops || num != 1, "Wrong cycle size"); + } + + return; +} + + +int main() { + + for (int i = 0; i < lgfn; ++i) { + SmartGraph graph; + SmartGraph::EdgeMap weight(graph); + + istringstream lgfs(lgf[i]); + graphReader(graph, lgfs). + edgeMap("weight", weight).run(); + + bool perfect_with_loops; + { + MaxFractionalMatching mfm(graph, true); + mfm.run(); + checkFractionalMatching(graph, mfm, true); + perfect_with_loops = mfm.matchingSize() == countNodes(graph); + } + + bool perfect_without_loops; + { + MaxFractionalMatching mfm(graph, false); + mfm.run(); + checkFractionalMatching(graph, mfm, false); + perfect_without_loops = mfm.matchingSize() == countNodes(graph); + } + + { + MaxFractionalMatching mfm(graph, true); + bool result = mfm.runPerfect(); + checkPerfectFractionalMatching(graph, mfm, result, true); + check(result == perfect_with_loops, "Wrong perfect matching"); + } + + { + MaxFractionalMatching mfm(graph, false); + bool result = mfm.runPerfect(); + checkPerfectFractionalMatching(graph, mfm, result, false); + check(result == perfect_without_loops, "Wrong perfect matching"); + } + + { + MaxWeightedFractionalMatching mwfm(graph, weight, true); + mwfm.run(); + checkWeightedFractionalMatching(graph, weight, mwfm, true); + } + + { + MaxWeightedFractionalMatching mwfm(graph, weight, false); + mwfm.run(); + checkWeightedFractionalMatching(graph, weight, mwfm, false); + } + + { + MaxWeightedPerfectFractionalMatching mwpfm(graph, weight, + true); + bool perfect = mwpfm.run(); + check(perfect == (mwpfm.matchingSize() == countNodes(graph)), + "Perfect matching found"); + check(perfect == perfect_with_loops, "Wrong perfect matching"); + + if (perfect) { + checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, true); + } + } + + { + MaxWeightedPerfectFractionalMatching mwpfm(graph, weight, + false); + bool perfect = mwpfm.run(); + check(perfect == (mwpfm.matchingSize() == countNodes(graph)), + "Perfect matching found"); + check(perfect == perfect_without_loops, "Wrong perfect matching"); + + if (perfect) { + checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, false); + } + } + + } + + return 0; +}