/* -*- 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