diff --git a/doc/groups.dox b/doc/groups.dox --- a/doc/groups.dox +++ b/doc/groups.dox @@ -386,7 +386,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. @@ -522,6 +522,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 matching.png \image latex matching.eps "Min Cost Perfect Matching" width=\textwidth diff --git a/lemon/Makefile.am b/lemon/Makefile.am --- a/lemon/Makefile.am +++ b/lemon/Makefile.am @@ -84,6 +84,7 @@ lemon/error.h \ lemon/euler.h \ lemon/fib_heap.h \ + lemon/fractional_matching.h \ lemon/full_graph.h \ lemon/glpk.h \ lemon/gomory_hu.h \ diff --git a/lemon/fractional_matching.h b/lemon/fractional_matching.h new file mode 100644 --- /dev/null +++ b/lemon/fractional_matching.h @@ -0,0 +1,2130 @@ +/* -*- 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 MaxFractionalMatching::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 uses 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. + /// + /// The primal solution is multiplied by + /// \ref MaxWeightedFractionalMatching::primalScale "2". + /// If the value type is integer, then the dual + /// solution is scaled by + /// \ref MaxWeightedFractionalMatching::dualScale "4". + /// + /// \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. + static const int primalScale = 2; + + /// \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 %MaxWeightedFractionalMatching 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. + 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 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 uses 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. + /// + /// The primal solution is multiplied by + /// \ref MaxWeightedPerfectFractionalMatching::primalScale "2". + /// If the value type is integer, then the dual + /// solution is scaled by + /// \ref MaxWeightedPerfectFractionalMatching::dualScale "4". + /// + /// \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. + static const int primalScale = 2; + + /// \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 %MaxWeightedPerfectFractionalMatching + /// 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. + 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 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 --git a/lemon/matching.h b/lemon/matching.h --- a/lemon/matching.h +++ b/lemon/matching.h @@ -16,8 +16,8 @@ * */ -#ifndef LEMON_MAX_MATCHING_H -#define LEMON_MAX_MATCHING_H +#ifndef LEMON_MATCHING_H +#define LEMON_MATCHING_H #include #include @@ -28,6 +28,7 @@ #include #include #include +#include ///\ingroup matching ///\file @@ -41,7 +42,7 @@ /// /// This class implements Edmonds' alternating forest matching algorithm /// for finding a maximum cardinality matching in a general undirected graph. - /// It can be started from an arbitrary initial matching + /// It can be started from an arbitrary initial matching /// (the default is the empty one). /// /// The dual solution of the problem is a map of the nodes to @@ -69,11 +70,11 @@ ///\brief Status constants for Gallai-Edmonds decomposition. /// - ///These constants are used for indicating the Gallai-Edmonds + ///These constants are used for indicating the Gallai-Edmonds ///decomposition of a graph. The nodes with status \c EVEN (or \c D) ///induce a subgraph with factor-critical components, the nodes with ///status \c ODD (or \c A) form the canonical barrier, and the nodes - ///with status \c MATCHED (or \c C) induce a subgraph having a + ///with status \c MATCHED (or \c C) induce a subgraph having a ///perfect matching. enum Status { EVEN = 1, ///< = 1. (\c D is an alias for \c EVEN.) @@ -512,7 +513,7 @@ } } - /// \brief Start Edmonds' algorithm with a heuristic improvement + /// \brief Start Edmonds' algorithm with a heuristic improvement /// for dense graphs /// /// This function runs Edmonds' algorithm with a heuristic of postponing @@ -534,8 +535,8 @@ /// \brief Run Edmonds' algorithm /// - /// This function runs Edmonds' algorithm. An additional heuristic of - /// postponing shrinks is used for relatively dense graphs + /// This function runs Edmonds' algorithm. An additional heuristic of + /// postponing shrinks is used for relatively dense graphs /// (for which m>=2*n holds). void run() { if (countEdges(_graph) < 2 * countNodes(_graph)) { @@ -556,7 +557,7 @@ /// \brief Return the size (cardinality) of the matching. /// - /// This function returns the size (cardinality) of the current matching. + /// This function returns the size (cardinality) of the current matching. /// After run() it returns the size of the maximum matching in the graph. int matchingSize() const { int size = 0; @@ -570,7 +571,7 @@ /// \brief Return \c true if the given edge is in the matching. /// - /// This function returns \c true if the given edge is in the current + /// This function returns \c true if the given edge is in the current /// matching. bool matching(const Edge& edge) const { return edge == (*_matching)[_graph.u(edge)]; @@ -579,7 +580,7 @@ /// \brief Return the matching arc (or edge) incident to the given node. /// /// This function returns the matching arc (or edge) incident to the - /// given node in the current matching or \c INVALID if the node is + /// given node in the current matching or \c INVALID if the node is /// not covered by the matching. Arc matching(const Node& n) const { return (*_matching)[n]; @@ -595,7 +596,7 @@ /// \brief Return the mate of the given node. /// - /// This function returns the mate of the given node in the current + /// This function returns the mate of the given node in the current /// matching or \c INVALID if the node is not covered by the matching. Node mate(const Node& n) const { return (*_matching)[n] != INVALID ? @@ -605,7 +606,7 @@ /// @} /// \name Dual Solution - /// Functions to get the dual solution, i.e. the Gallai-Edmonds + /// Functions to get the dual solution, i.e. the Gallai-Edmonds /// decomposition. /// @{ @@ -648,8 +649,8 @@ /// on extensive use of priority queues and provides /// \f$O(nm\log n)\f$ time complexity. /// - /// The maximum weighted matching problem is to find a subset of the - /// edges in an undirected graph with maximum overall weight for which + /// The maximum weighted matching problem is to find a subset of the + /// edges in an undirected graph with maximum overall weight for which /// each node has at most one incident edge. /// 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] @@ -673,16 +674,16 @@ /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} \frac{\vert B \vert - 1}{2}z_B\f] */ /// - /// The algorithm can be executed with the run() function. + /// 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 and the - /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, - /// which is able to iterate on the nodes of a blossom. + /// can be obtained using the query functions and the + /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, + /// which is able to iterate on the nodes of a blossom. /// If the value type is integer, then the dual solution is multiplied /// by \ref MaxWeightedMatching::dualScale "4". /// /// \tparam GR The undirected graph type the algorithm runs on. - /// \tparam WM The type edge weight map. The default type is + /// \tparam WM The type edge weight map. The default type is /// \ref concepts::Graph::EdgeMap "GR::EdgeMap". #ifdef DOXYGEN template @@ -745,7 +746,7 @@ typedef RangeMap IntIntMap; enum Status { - EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2 + EVEN = -1, MATCHED = 0, ODD = 1 }; typedef HeapUnionFind BlossomSet; @@ -797,6 +798,10 @@ BinHeap *_delta4; Value _delta_sum; + int _unmatched; + + typedef MaxWeightedFractionalMatching FractionalMatching; + FractionalMatching *_fractional; void createStructures() { _node_num = countNodes(_graph); @@ -844,9 +849,6 @@ } void destroyStructures() { - _node_num = countNodes(_graph); - _blossom_num = _node_num * 3 / 2; - if (_matching) { delete _matching; } @@ -922,10 +924,6 @@ if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { _delta3->push(e, rw / 2); } - } else if ((*_blossom_data)[vb].status == UNMATCHED) { - if (_delta3->state(e) != _delta3->IN_HEAP) { - _delta3->push(e, rw); - } } else { typename std::map::iterator it = (*_node_data)[vi].heap_index.find(tree); @@ -949,202 +947,6 @@ _delta2->push(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset){ - _delta2->decrease(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } - } - } - } - } - } - (*_blossom_data)[blossom].offset = 0; - } - - void matchedToOdd(int blossom) { - if (_delta2->state(blossom) == _delta2->IN_HEAP) { - _delta2->erase(blossom); - } - (*_blossom_data)[blossom].offset += _delta_sum; - if (!_blossom_set->trivial(blossom)) { - _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + - (*_blossom_data)[blossom].offset); - } - } - - void evenToMatched(int blossom, int tree) { - if (!_blossom_set->trivial(blossom)) { - (*_blossom_data)[blossom].pot += 2 * _delta_sum; - } - - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); - n != INVALID; ++n) { - int ni = (*_node_index)[n]; - (*_node_data)[ni].pot -= _delta_sum; - - _delta1->erase(n); - - for (InArcIt e(_graph, n); e != INVALID; ++e) { - Node v = _graph.source(e); - int vb = _blossom_set->find(v); - int vi = (*_node_index)[v]; - - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - - dualScale * _weight[e]; - - if (vb == blossom) { - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } - } else if ((*_blossom_data)[vb].status == EVEN) { - - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } - - int vt = _tree_set->find(vb); - - if (vt != tree) { - - Arc r = _graph.oppositeArc(e); - - typename std::map::iterator it = - (*_node_data)[ni].heap_index.find(vt); - - if (it != (*_node_data)[ni].heap_index.end()) { - if ((*_node_data)[ni].heap[it->second] > rw) { - (*_node_data)[ni].heap.replace(it->second, r); - (*_node_data)[ni].heap.decrease(r, rw); - it->second = r; - } - } else { - (*_node_data)[ni].heap.push(r, rw); - (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); - } - - if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { - _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); - - if (_delta2->state(blossom) != _delta2->IN_HEAP) { - _delta2->push(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); - } else if ((*_delta2)[blossom] > - _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset){ - _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); - } - } - } - - } else if ((*_blossom_data)[vb].status == UNMATCHED) { - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } - } else { - - typename std::map::iterator it = - (*_node_data)[vi].heap_index.find(tree); - - if (it != (*_node_data)[vi].heap_index.end()) { - (*_node_data)[vi].heap.erase(it->second); - (*_node_data)[vi].heap_index.erase(it); - if ((*_node_data)[vi].heap.empty()) { - _blossom_set->increase(v, std::numeric_limits::max()); - } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { - _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); - } - - if ((*_blossom_data)[vb].status == MATCHED) { - if (_blossom_set->classPrio(vb) == - std::numeric_limits::max()) { - _delta2->erase(vb); - } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset) { - _delta2->increase(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } - } - } - } - } - } - } - - void oddToMatched(int blossom) { - (*_blossom_data)[blossom].offset -= _delta_sum; - - if (_blossom_set->classPrio(blossom) != - std::numeric_limits::max()) { - _delta2->push(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); - } - - if (!_blossom_set->trivial(blossom)) { - _delta4->erase(blossom); - } - } - - void oddToEven(int blossom, int tree) { - if (!_blossom_set->trivial(blossom)) { - _delta4->erase(blossom); - (*_blossom_data)[blossom].pot -= - 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); - } - - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); - n != INVALID; ++n) { - int ni = (*_node_index)[n]; - - _blossom_set->increase(n, std::numeric_limits::max()); - - (*_node_data)[ni].heap.clear(); - (*_node_data)[ni].heap_index.clear(); - (*_node_data)[ni].pot += - 2 * _delta_sum - (*_blossom_data)[blossom].offset; - - _delta1->push(n, (*_node_data)[ni].pot); - - for (InArcIt e(_graph, n); e != INVALID; ++e) { - Node v = _graph.source(e); - int vb = _blossom_set->find(v); - int vi = (*_node_index)[v]; - - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - - dualScale * _weight[e]; - - if ((*_blossom_data)[vb].status == EVEN) { - if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { - _delta3->push(e, rw / 2); - } - } else if ((*_blossom_data)[vb].status == UNMATCHED) { - if (_delta3->state(e) != _delta3->IN_HEAP) { - _delta3->push(e, rw); - } - } else { - - typename std::map::iterator it = - (*_node_data)[vi].heap_index.find(tree); - - if (it != (*_node_data)[vi].heap_index.end()) { - if ((*_node_data)[vi].heap[it->second] > rw) { - (*_node_data)[vi].heap.replace(it->second, e); - (*_node_data)[vi].heap.decrease(e, rw); - it->second = e; - } - } else { - (*_node_data)[vi].heap.push(e, rw); - (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); - } - - if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { - _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); - - if ((*_blossom_data)[vb].status == MATCHED) { - if (_delta2->state(vb) != _delta2->IN_HEAP) { - _delta2->push(vb, _blossom_set->classPrio(vb) - - (*_blossom_data)[vb].offset); - } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset) { _delta2->decrease(vb, _blossom_set->classPrio(vb) - (*_blossom_data)[vb].offset); @@ -1157,43 +959,145 @@ (*_blossom_data)[blossom].offset = 0; } - - void matchedToUnmatched(int blossom) { + void matchedToOdd(int blossom) { if (_delta2->state(blossom) == _delta2->IN_HEAP) { _delta2->erase(blossom); } + (*_blossom_data)[blossom].offset += _delta_sum; + if (!_blossom_set->trivial(blossom)) { + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + + (*_blossom_data)[blossom].offset); + } + } + + void evenToMatched(int blossom, int tree) { + if (!_blossom_set->trivial(blossom)) { + (*_blossom_data)[blossom].pot += 2 * _delta_sum; + } for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); n != INVALID; ++n) { int ni = (*_node_index)[n]; - - _blossom_set->increase(n, std::numeric_limits::max()); - - (*_node_data)[ni].heap.clear(); - (*_node_data)[ni].heap_index.clear(); - - for (OutArcIt e(_graph, n); e != INVALID; ++e) { - Node v = _graph.target(e); + (*_node_data)[ni].pot -= _delta_sum; + + _delta1->erase(n); + + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Node v = _graph.source(e); int vb = _blossom_set->find(v); int vi = (*_node_index)[v]; Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - dualScale * _weight[e]; - if ((*_blossom_data)[vb].status == EVEN) { - if (_delta3->state(e) != _delta3->IN_HEAP) { - _delta3->push(e, rw); + if (vb == blossom) { + if (_delta3->state(e) == _delta3->IN_HEAP) { + _delta3->erase(e); + } + } else if ((*_blossom_data)[vb].status == EVEN) { + + if (_delta3->state(e) == _delta3->IN_HEAP) { + _delta3->erase(e); + } + + int vt = _tree_set->find(vb); + + if (vt != tree) { + + Arc r = _graph.oppositeArc(e); + + typename std::map::iterator it = + (*_node_data)[ni].heap_index.find(vt); + + if (it != (*_node_data)[ni].heap_index.end()) { + if ((*_node_data)[ni].heap[it->second] > rw) { + (*_node_data)[ni].heap.replace(it->second, r); + (*_node_data)[ni].heap.decrease(r, rw); + it->second = r; + } + } else { + (*_node_data)[ni].heap.push(r, rw); + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); + } + + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); + + if (_delta2->state(blossom) != _delta2->IN_HEAP) { + _delta2->push(blossom, _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset); + } else if ((*_delta2)[blossom] > + _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset){ + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset); + } + } + } + } else { + + typename std::map::iterator it = + (*_node_data)[vi].heap_index.find(tree); + + if (it != (*_node_data)[vi].heap_index.end()) { + (*_node_data)[vi].heap.erase(it->second); + (*_node_data)[vi].heap_index.erase(it); + if ((*_node_data)[vi].heap.empty()) { + _blossom_set->increase(v, std::numeric_limits::max()); + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { + _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); + } + + if ((*_blossom_data)[vb].status == MATCHED) { + if (_blossom_set->classPrio(vb) == + std::numeric_limits::max()) { + _delta2->erase(vb); + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset) { + _delta2->increase(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } + } } } } } } - void unmatchedToMatched(int blossom) { + void oddToMatched(int blossom) { + (*_blossom_data)[blossom].offset -= _delta_sum; + + if (_blossom_set->classPrio(blossom) != + std::numeric_limits::max()) { + _delta2->push(blossom, _blossom_set->classPrio(blossom) - + (*_blossom_data)[blossom].offset); + } + + if (!_blossom_set->trivial(blossom)) { + _delta4->erase(blossom); + } + } + + void oddToEven(int blossom, int tree) { + if (!_blossom_set->trivial(blossom)) { + _delta4->erase(blossom); + (*_blossom_data)[blossom].pot -= + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); + } + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); n != INVALID; ++n) { int ni = (*_node_index)[n]; + _blossom_set->increase(n, std::numeric_limits::max()); + + (*_node_data)[ni].heap.clear(); + (*_node_data)[ni].heap_index.clear(); + (*_node_data)[ni].pot += + 2 * _delta_sum - (*_blossom_data)[blossom].offset; + + _delta1->push(n, (*_node_data)[ni].pot); + for (InArcIt e(_graph, n); e != INVALID; ++e) { Node v = _graph.source(e); int vb = _blossom_set->find(v); @@ -1202,54 +1106,44 @@ Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - dualScale * _weight[e]; - if (vb == blossom) { - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); + if ((*_blossom_data)[vb].status == EVEN) { + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { + _delta3->push(e, rw / 2); } - } else if ((*_blossom_data)[vb].status == EVEN) { - - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } - - int vt = _tree_set->find(vb); - - Arc r = _graph.oppositeArc(e); + } else { typename std::map::iterator it = - (*_node_data)[ni].heap_index.find(vt); - - if (it != (*_node_data)[ni].heap_index.end()) { - if ((*_node_data)[ni].heap[it->second] > rw) { - (*_node_data)[ni].heap.replace(it->second, r); - (*_node_data)[ni].heap.decrease(r, rw); - it->second = r; + (*_node_data)[vi].heap_index.find(tree); + + if (it != (*_node_data)[vi].heap_index.end()) { + if ((*_node_data)[vi].heap[it->second] > rw) { + (*_node_data)[vi].heap.replace(it->second, e); + (*_node_data)[vi].heap.decrease(e, rw); + it->second = e; } } else { - (*_node_data)[ni].heap.push(r, rw); - (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); + (*_node_data)[vi].heap.push(e, rw); + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); } - if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { - _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); - - if (_delta2->state(blossom) != _delta2->IN_HEAP) { - _delta2->push(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); - } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)- - (*_blossom_data)[blossom].offset){ - _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - - (*_blossom_data)[blossom].offset); + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); + + if ((*_blossom_data)[vb].status == MATCHED) { + if (_delta2->state(vb) != _delta2->IN_HEAP) { + _delta2->push(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset) { + _delta2->decrease(vb, _blossom_set->classPrio(vb) - + (*_blossom_data)[vb].offset); + } } } - - } else if ((*_blossom_data)[vb].status == UNMATCHED) { - if (_delta3->state(e) == _delta3->IN_HEAP) { - _delta3->erase(e); - } } } } + (*_blossom_data)[blossom].offset = 0; } void alternatePath(int even, int tree) { @@ -1294,39 +1188,42 @@ alternatePath(blossom, tree); destroyTree(tree); - (*_blossom_data)[blossom].status = UNMATCHED; (*_blossom_data)[blossom].base = node; - matchedToUnmatched(blossom); + (*_blossom_data)[blossom].next = INVALID; } - void augmentOnEdge(const Edge& edge) { int left = _blossom_set->find(_graph.u(edge)); int right = _blossom_set->find(_graph.v(edge)); - if ((*_blossom_data)[left].status == EVEN) { - int left_tree = _tree_set->find(left); - alternatePath(left, left_tree); - destroyTree(left_tree); - } else { - (*_blossom_data)[left].status = MATCHED; - unmatchedToMatched(left); - } - - if ((*_blossom_data)[right].status == EVEN) { - int right_tree = _tree_set->find(right); - alternatePath(right, right_tree); - destroyTree(right_tree); - } else { - (*_blossom_data)[right].status = MATCHED; - unmatchedToMatched(right); - } + int left_tree = _tree_set->find(left); + alternatePath(left, left_tree); + destroyTree(left_tree); + + int right_tree = _tree_set->find(right); + alternatePath(right, right_tree); + destroyTree(right_tree); (*_blossom_data)[left].next = _graph.direct(edge, true); (*_blossom_data)[right].next = _graph.direct(edge, false); } + void augmentOnArc(const Arc& arc) { + + int left = _blossom_set->find(_graph.source(arc)); + int right = _blossom_set->find(_graph.target(arc)); + + (*_blossom_data)[left].status = MATCHED; + + int right_tree = _tree_set->find(right); + alternatePath(right, right_tree); + destroyTree(right_tree); + + (*_blossom_data)[left].next = arc; + (*_blossom_data)[right].next = _graph.oppositeArc(arc); + } + void extendOnArc(const Arc& arc) { int base = _blossom_set->find(_graph.target(arc)); int tree = _tree_set->find(base); @@ -1529,7 +1426,7 @@ _tree_set->insert(sb, tree); (*_blossom_data)[sb].pred = pred; (*_blossom_data)[sb].next = - _graph.oppositeArc((*_blossom_data)[tb].next); + _graph.oppositeArc((*_blossom_data)[tb].next); pred = (*_blossom_data)[ub].next; @@ -1629,7 +1526,7 @@ } for (int i = 0; i < int(blossoms.size()); ++i) { - if ((*_blossom_data)[blossoms[i]].status == MATCHED) { + if ((*_blossom_data)[blossoms[i]].next != INVALID) { Value offset = (*_blossom_data)[blossoms[i]].offset; (*_blossom_data)[blossoms[i]].pot += 2 * offset; @@ -1667,10 +1564,16 @@ _delta3_index(0), _delta3(0), _delta4_index(0), _delta4(0), - _delta_sum() {} + _delta_sum(), _unmatched(0), + + _fractional(0) + {} ~MaxWeightedMatching() { destroyStructures(); + if (_fractional) { + delete _fractional; + } } /// \name Execution Control @@ -1699,6 +1602,8 @@ (*_delta4_index)[i] = _delta4->PRE_HEAP; } + _unmatched = _node_num; + int index = 0; for (NodeIt n(_graph); n != INVALID; ++n) { Value max = 0; @@ -1733,18 +1638,155 @@ } } + /// \brief Initialize the algorithm with fractional matching + /// + /// This function initializes the algorithm with a fractional + /// matching. This initialization is also called jumpstart heuristic. + void fractionalInit() { + createStructures(); + + if (_fractional == 0) { + _fractional = new FractionalMatching(_graph, _weight, false); + } + _fractional->run(); + + for (ArcIt e(_graph); e != INVALID; ++e) { + (*_node_heap_index)[e] = BinHeap::PRE_HEAP; + } + for (NodeIt n(_graph); n != INVALID; ++n) { + (*_delta1_index)[n] = _delta1->PRE_HEAP; + } + for (EdgeIt e(_graph); e != INVALID; ++e) { + (*_delta3_index)[e] = _delta3->PRE_HEAP; + } + for (int i = 0; i < _blossom_num; ++i) { + (*_delta2_index)[i] = _delta2->PRE_HEAP; + (*_delta4_index)[i] = _delta4->PRE_HEAP; + } + + _unmatched = 0; + + int index = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + Value pot = _fractional->nodeValue(n); + (*_node_index)[n] = index; + (*_node_data)[index].pot = pot; + int blossom = + _blossom_set->insert(n, std::numeric_limits::max()); + + (*_blossom_data)[blossom].status = MATCHED; + (*_blossom_data)[blossom].pred = INVALID; + (*_blossom_data)[blossom].next = _fractional->matching(n); + if (_fractional->matching(n) == INVALID) { + (*_blossom_data)[blossom].base = n; + } + (*_blossom_data)[blossom].pot = 0; + (*_blossom_data)[blossom].offset = 0; + ++index; + } + + typename Graph::template NodeMap processed(_graph, false); + for (NodeIt n(_graph); n != INVALID; ++n) { + if (processed[n]) continue; + processed[n] = true; + if (_fractional->matching(n) == INVALID) continue; + int num = 1; + Node v = _graph.target(_fractional->matching(n)); + while (n != v) { + processed[v] = true; + v = _graph.target(_fractional->matching(v)); + ++num; + } + + if (num % 2 == 1) { + std::vector subblossoms(num); + + subblossoms[--num] = _blossom_set->find(n); + _delta1->push(n, _fractional->nodeValue(n)); + v = _graph.target(_fractional->matching(n)); + while (n != v) { + subblossoms[--num] = _blossom_set->find(v); + _delta1->push(v, _fractional->nodeValue(v)); + v = _graph.target(_fractional->matching(v)); + } + + int surface = + _blossom_set->join(subblossoms.begin(), subblossoms.end()); + (*_blossom_data)[surface].status = EVEN; + (*_blossom_data)[surface].pred = INVALID; + (*_blossom_data)[surface].next = INVALID; + (*_blossom_data)[surface].pot = 0; + (*_blossom_data)[surface].offset = 0; + + _tree_set->insert(surface); + ++_unmatched; + } + } + + for (EdgeIt e(_graph); e != INVALID; ++e) { + int si = (*_node_index)[_graph.u(e)]; + int sb = _blossom_set->find(_graph.u(e)); + int ti = (*_node_index)[_graph.v(e)]; + int tb = _blossom_set->find(_graph.v(e)); + if ((*_blossom_data)[sb].status == EVEN && + (*_blossom_data)[tb].status == EVEN && sb != tb) { + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - + dualScale * _weight[e]) / 2); + } + } + + for (NodeIt n(_graph); n != INVALID; ++n) { + int nb = _blossom_set->find(n); + if ((*_blossom_data)[nb].status != MATCHED) continue; + int ni = (*_node_index)[n]; + + for (OutArcIt e(_graph, n); e != INVALID; ++e) { + Node v = _graph.target(e); + int vb = _blossom_set->find(v); + int vi = (*_node_index)[v]; + + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - + dualScale * _weight[e]; + + if ((*_blossom_data)[vb].status == EVEN) { + + int vt = _tree_set->find(vb); + + typename std::map::iterator it = + (*_node_data)[ni].heap_index.find(vt); + + if (it != (*_node_data)[ni].heap_index.end()) { + if ((*_node_data)[ni].heap[it->second] > rw) { + (*_node_data)[ni].heap.replace(it->second, e); + (*_node_data)[ni].heap.decrease(e, rw); + it->second = e; + } + } else { + (*_node_data)[ni].heap.push(e, rw); + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); + } + } + } + + if (!(*_node_data)[ni].heap.empty()) { + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); + _delta2->push(nb, _blossom_set->classPrio(nb)); + } + } + } + /// \brief Start the algorithm /// /// This function starts the algorithm. /// - /// \pre \ref init() must be called before using this function. + /// \pre \ref init() or \ref fractionalInit() must be called + /// before using this function. void start() { enum OpType { D1, D2, D3, D4 }; - int unmatched = _node_num; - while (unmatched > 0) { + while (_unmatched > 0) { Value d1 = !_delta1->empty() ? _delta1->prio() : std::numeric_limits::max(); @@ -1757,26 +1799,30 @@ Value d4 = !_delta4->empty() ? _delta4->prio() : std::numeric_limits::max(); - _delta_sum = d1; OpType ot = D1; + _delta_sum = d3; OpType ot = D3; + if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; } if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } - if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } - switch (ot) { case D1: { Node n = _delta1->top(); unmatchNode(n); - --unmatched; + --_unmatched; } break; case D2: { int blossom = _delta2->top(); Node n = _blossom_set->classTop(blossom); - Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); - extendOnArc(e); + Arc a = (*_node_data)[(*_node_index)[n]].heap.top(); + if ((*_blossom_data)[blossom].next == INVALID) { + augmentOnArc(a); + --_unmatched; + } else { + extendOnArc(a); + } } break; case D3: @@ -1789,26 +1835,14 @@ if (left_blossom == right_blossom) { _delta3->pop(); } else { - int left_tree; - if ((*_blossom_data)[left_blossom].status == EVEN) { - left_tree = _tree_set->find(left_blossom); - } else { - left_tree = -1; - ++unmatched; - } - int right_tree; - if ((*_blossom_data)[right_blossom].status == EVEN) { - right_tree = _tree_set->find(right_blossom); - } else { - right_tree = -1; - ++unmatched; - } + int left_tree = _tree_set->find(left_blossom); + int right_tree = _tree_set->find(right_blossom); if (left_tree == right_tree) { shrinkOnEdge(e, left_tree); } else { augmentOnEdge(e); - unmatched -= 2; + _unmatched -= 2; } } } break; @@ -1826,18 +1860,18 @@ /// /// \note mwm.run() is just a shortcut of the following code. /// \code - /// mwm.init(); + /// mwm.fractionalInit(); /// mwm.start(); /// \endcode void run() { - init(); + fractionalInit(); start(); } /// @} /// \name Primal Solution - /// Functions to get the primal solution, i.e. the maximum weighted + /// 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. @@ -1856,7 +1890,7 @@ sum += _weight[(*_matching)[n]]; } } - return sum /= 2; + return sum / 2; } /// \brief Return the size (cardinality) of the matching. @@ -1876,7 +1910,7 @@ /// \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 + /// This function returns \c true if the given edge is in the found /// matching. /// /// \pre Either run() or start() must be called before using this function. @@ -1887,7 +1921,7 @@ /// \brief Return the matching arc (or edge) incident to the given node. /// /// This function returns the matching arc (or edge) incident to the - /// given node in the found matching or \c INVALID if the node is + /// given node in the found matching or \c INVALID if the node is /// not covered by the matching. /// /// \pre Either run() or start() must be called before using this function. @@ -1905,7 +1939,7 @@ /// \brief Return the mate of the given node. /// - /// This function returns the mate of the given node in the found + /// This function returns the mate of the given node in the found /// matching or \c INVALID if the node is not covered by the matching. /// /// \pre Either run() or start() must be called before using this function. @@ -1925,8 +1959,8 @@ /// \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 + /// 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. @@ -1981,9 +2015,9 @@ /// \brief Iterator for obtaining the nodes of a blossom. /// - /// This class provides an iterator for obtaining the nodes of the + /// This class provides an iterator for obtaining the nodes of the /// given blossom. It lists a subset of the nodes. - /// Before using this iterator, you must allocate a + /// Before using this iterator, you must allocate a /// MaxWeightedMatching class and execute it. class BlossomIt { public: @@ -1992,8 +2026,8 @@ /// /// Constructor to get the nodes of the given variable. /// - /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or - /// \ref MaxWeightedMatching::start() "algorithm.start()" must be + /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or + /// \ref MaxWeightedMatching::start() "algorithm.start()" must be /// called before initializing this iterator. BlossomIt(const MaxWeightedMatching& algorithm, int variable) : _algorithm(&algorithm) @@ -2046,8 +2080,8 @@ /// is based on extensive use of priority queues and provides /// \f$O(nm\log n)\f$ time complexity. /// - /// The maximum weighted perfect matching problem is to find a subset of - /// the edges in an undirected graph with maximum overall weight for which + /// The maximum weighted perfect matching problem is to find a subset of + /// the edges in an undirected graph with maximum overall weight for which /// each node has exactly one incident edge. /// It can be formulated with the following linear program. /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f] @@ -2070,16 +2104,16 @@ /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} \frac{\vert B \vert - 1}{2}z_B\f] */ /// - /// The algorithm can be executed with the run() function. + /// 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 and the - /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, - /// which is able to iterate on the nodes of a blossom. + /// can be obtained using the query functions and the + /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, + /// which is able to iterate on the nodes of a blossom. /// If the value type is integer, then the dual solution is multiplied /// by \ref MaxWeightedMatching::dualScale "4". /// /// \tparam GR The undirected graph type the algorithm runs on. - /// \tparam WM The type edge weight map. The default type is + /// \tparam WM The type edge weight map. The default type is /// \ref concepts::Graph::EdgeMap "GR::EdgeMap". #ifdef DOXYGEN template @@ -2190,6 +2224,11 @@ BinHeap *_delta4; Value _delta_sum; + int _unmatched; + + typedef MaxWeightedPerfectFractionalMatching + FractionalMatching; + FractionalMatching *_fractional; void createStructures() { _node_num = countNodes(_graph); @@ -2233,9 +2272,6 @@ } void destroyStructures() { - _node_num = countNodes(_graph); - _blossom_num = _node_num * 3 / 2; - if (_matching) { delete _matching; } @@ -2908,10 +2944,16 @@ _delta3_index(0), _delta3(0), _delta4_index(0), _delta4(0), - _delta_sum() {} + _delta_sum(), _unmatched(0), + + _fractional(0) + {} ~MaxWeightedPerfectMatching() { destroyStructures(); + if (_fractional) { + delete _fractional; + } } /// \name Execution Control @@ -2937,6 +2979,8 @@ (*_delta4_index)[i] = _delta4->PRE_HEAP; } + _unmatched = _node_num; + int index = 0; for (NodeIt n(_graph); n != INVALID; ++n) { Value max = - std::numeric_limits::max(); @@ -2970,18 +3014,152 @@ } } + /// \brief Initialize the algorithm with fractional matching + /// + /// This function initializes the algorithm with a fractional + /// matching. This initialization is also called jumpstart heuristic. + void fractionalInit() { + createStructures(); + + if (_fractional == 0) { + _fractional = new FractionalMatching(_graph, _weight, false); + } + if (!_fractional->run()) { + _unmatched = -1; + return; + } + + for (ArcIt e(_graph); e != INVALID; ++e) { + (*_node_heap_index)[e] = BinHeap::PRE_HEAP; + } + for (EdgeIt e(_graph); e != INVALID; ++e) { + (*_delta3_index)[e] = _delta3->PRE_HEAP; + } + for (int i = 0; i < _blossom_num; ++i) { + (*_delta2_index)[i] = _delta2->PRE_HEAP; + (*_delta4_index)[i] = _delta4->PRE_HEAP; + } + + _unmatched = 0; + + int index = 0; + for (NodeIt n(_graph); n != INVALID; ++n) { + Value pot = _fractional->nodeValue(n); + (*_node_index)[n] = index; + (*_node_data)[index].pot = pot; + int blossom = + _blossom_set->insert(n, std::numeric_limits::max()); + + (*_blossom_data)[blossom].status = MATCHED; + (*_blossom_data)[blossom].pred = INVALID; + (*_blossom_data)[blossom].next = _fractional->matching(n); + (*_blossom_data)[blossom].pot = 0; + (*_blossom_data)[blossom].offset = 0; + ++index; + } + + typename Graph::template NodeMap processed(_graph, false); + for (NodeIt n(_graph); n != INVALID; ++n) { + if (processed[n]) continue; + processed[n] = true; + if (_fractional->matching(n) == INVALID) continue; + int num = 1; + Node v = _graph.target(_fractional->matching(n)); + while (n != v) { + processed[v] = true; + v = _graph.target(_fractional->matching(v)); + ++num; + } + + if (num % 2 == 1) { + std::vector subblossoms(num); + + subblossoms[--num] = _blossom_set->find(n); + v = _graph.target(_fractional->matching(n)); + while (n != v) { + subblossoms[--num] = _blossom_set->find(v); + v = _graph.target(_fractional->matching(v)); + } + + int surface = + _blossom_set->join(subblossoms.begin(), subblossoms.end()); + (*_blossom_data)[surface].status = EVEN; + (*_blossom_data)[surface].pred = INVALID; + (*_blossom_data)[surface].next = INVALID; + (*_blossom_data)[surface].pot = 0; + (*_blossom_data)[surface].offset = 0; + + _tree_set->insert(surface); + ++_unmatched; + } + } + + for (EdgeIt e(_graph); e != INVALID; ++e) { + int si = (*_node_index)[_graph.u(e)]; + int sb = _blossom_set->find(_graph.u(e)); + int ti = (*_node_index)[_graph.v(e)]; + int tb = _blossom_set->find(_graph.v(e)); + if ((*_blossom_data)[sb].status == EVEN && + (*_blossom_data)[tb].status == EVEN && sb != tb) { + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - + dualScale * _weight[e]) / 2); + } + } + + for (NodeIt n(_graph); n != INVALID; ++n) { + int nb = _blossom_set->find(n); + if ((*_blossom_data)[nb].status != MATCHED) continue; + int ni = (*_node_index)[n]; + + for (OutArcIt e(_graph, n); e != INVALID; ++e) { + Node v = _graph.target(e); + int vb = _blossom_set->find(v); + int vi = (*_node_index)[v]; + + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - + dualScale * _weight[e]; + + if ((*_blossom_data)[vb].status == EVEN) { + + int vt = _tree_set->find(vb); + + typename std::map::iterator it = + (*_node_data)[ni].heap_index.find(vt); + + if (it != (*_node_data)[ni].heap_index.end()) { + if ((*_node_data)[ni].heap[it->second] > rw) { + (*_node_data)[ni].heap.replace(it->second, e); + (*_node_data)[ni].heap.decrease(e, rw); + it->second = e; + } + } else { + (*_node_data)[ni].heap.push(e, rw); + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); + } + } + } + + if (!(*_node_data)[ni].heap.empty()) { + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); + _delta2->push(nb, _blossom_set->classPrio(nb)); + } + } + } + /// \brief Start the algorithm /// /// This function starts the algorithm. /// - /// \pre \ref init() must be called before using this function. + /// \pre \ref init() or \ref fractionalInit() must be called before + /// using this function. bool start() { enum OpType { D2, D3, D4 }; - int unmatched = _node_num; - while (unmatched > 0) { + if (_unmatched == -1) return false; + + while (_unmatched > 0) { Value d2 = !_delta2->empty() ? _delta2->prio() : std::numeric_limits::max(); @@ -2991,8 +3169,8 @@ Value d4 = !_delta4->empty() ? _delta4->prio() : std::numeric_limits::max(); - _delta_sum = d2; OpType ot = D2; - if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } + _delta_sum = d3; OpType ot = D3; + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } if (_delta_sum == std::numeric_limits::max()) { @@ -3025,7 +3203,7 @@ shrinkOnEdge(e, left_tree); } else { augmentOnEdge(e); - unmatched -= 2; + _unmatched -= 2; } } } break; @@ -3044,18 +3222,18 @@ /// /// \note mwpm.run() is just a shortcut of the following code. /// \code - /// mwpm.init(); + /// mwpm.fractionalInit(); /// mwpm.start(); /// \endcode bool run() { - init(); + fractionalInit(); return start(); } /// @} /// \name Primal Solution - /// Functions to get the primal solution, i.e. the maximum weighted + /// Functions to get the primal solution, i.e. the maximum weighted /// perfect matching.\n /// Either \ref run() or \ref start() function should be called before /// using them. @@ -3074,12 +3252,12 @@ sum += _weight[(*_matching)[n]]; } } - return sum /= 2; + return sum / 2; } /// \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 + /// This function returns \c true if the given edge is in the found /// matching. /// /// \pre Either run() or start() must be called before using this function. @@ -3090,7 +3268,7 @@ /// \brief Return the matching arc (or edge) incident to the given node. /// /// This function returns the matching arc (or edge) incident to the - /// given node in the found matching or \c INVALID if the node is + /// given node in the found matching or \c INVALID if the node is /// not covered by the matching. /// /// \pre Either run() or start() must be called before using this function. @@ -3108,7 +3286,7 @@ /// \brief Return the mate of the given node. /// - /// This function returns the mate of the given node in the found + /// This function returns the mate of the given node in the found /// matching or \c INVALID if the node is not covered by the matching. /// /// \pre Either run() or start() must be called before using this function. @@ -3127,8 +3305,8 @@ /// \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 + /// 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. @@ -3183,9 +3361,9 @@ /// \brief Iterator for obtaining the nodes of a blossom. /// - /// This class provides an iterator for obtaining the nodes of the + /// This class provides an iterator for obtaining the nodes of the /// given blossom. It lists a subset of the nodes. - /// Before using this iterator, you must allocate a + /// Before using this iterator, you must allocate a /// MaxWeightedPerfectMatching class and execute it. class BlossomIt { public: @@ -3194,8 +3372,8 @@ /// /// Constructor to get the nodes of the given variable. /// - /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" - /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" + /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" + /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" /// must be called before initializing this iterator. BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) : _algorithm(&algorithm) @@ -3241,4 +3419,4 @@ } //END OF NAMESPACE LEMON -#endif //LEMON_MAX_MATCHING_H +#endif //LEMON_MATCHING_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -21,6 +21,7 @@ edge_set_test error_test euler_test + fractional_matching_test gomory_hu_test graph_copy_test graph_test diff --git a/test/Makefile.am b/test/Makefile.am --- a/test/Makefile.am +++ b/test/Makefile.am @@ -23,6 +23,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 \ @@ -71,6 +72,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 --git a/test/fractional_matching_test.cc b/test/fractional_matching_test.cc new file mode 100644 --- /dev/null +++ b/test/fractional_matching_test.cc @@ -0,0 +1,525 @@ +/* -*- 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"); + + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) { + check((e == mfm.matching(graph.u(e)) ? 1 : 0) + + (e == mfm.matching(graph.v(e)) ? 1 : 0) == + mfm.matching(e), "Invalid matching"); + } + + 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"); + } + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) { + check((e == mfm.matching(graph.u(e)) ? 1 : 0) + + (e == mfm.matching(graph.v(e)) ? 1 : 0) == + mfm.matching(e), "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"); + } + } + + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) { + check((e == mwfm.matching(graph.u(e)) ? 1 : 0) + + (e == mwfm.matching(graph.v(e)) ? 1 : 0) == + mwfm.matching(e), "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)); + } + + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) { + check((e == mwpfm.matching(graph.u(e)) ? 1 : 0) + + (e == mwpfm.matching(graph.v(e)) ? 1 : 0) == + mwpfm.matching(e), "Invalid matching"); + } + + 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; +} diff --git a/test/matching_test.cc b/test/matching_test.cc --- a/test/matching_test.cc +++ b/test/matching_test.cc @@ -401,22 +401,46 @@ graphReader(graph, lgfs). edgeMap("weight", weight).run(); - MaxMatching mm(graph); - mm.run(); - checkMatching(graph, mm); + bool perfect; + { + MaxMatching mm(graph); + mm.run(); + checkMatching(graph, mm); + perfect = 2 * mm.matchingSize() == countNodes(graph); + } - MaxWeightedMatching mwm(graph, weight); - mwm.run(); - checkWeightedMatching(graph, weight, mwm); + { + MaxWeightedMatching mwm(graph, weight); + mwm.run(); + checkWeightedMatching(graph, weight, mwm); + } - MaxWeightedPerfectMatching mwpm(graph, weight); - bool perfect = mwpm.run(); + { + MaxWeightedMatching mwm(graph, weight); + mwm.init(); + mwm.start(); + checkWeightedMatching(graph, weight, mwm); + } - check(perfect == (mm.matchingSize() * 2 == countNodes(graph)), - "Perfect matching found"); + { + MaxWeightedPerfectMatching mwpm(graph, weight); + bool result = mwpm.run(); + + check(result == perfect, "Perfect matching found"); + if (perfect) { + checkWeightedPerfectMatching(graph, weight, mwpm); + } + } - if (perfect) { - checkWeightedPerfectMatching(graph, weight, mwpm); + { + MaxWeightedPerfectMatching mwpm(graph, weight); + mwpm.init(); + bool result = mwpm.start(); + + check(result == perfect, "Perfect matching found"); + if (perfect) { + checkWeightedPerfectMatching(graph, weight, mwpm); + } } }