/* -*- C++ -*- * * This file is a part of LEMON, a generic C++ optimization library * * Copyright (C) 2003-2007 * 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_BIPARTITE_MATCHING #define LEMON_BIPARTITE_MATCHING #include #include #include #include ///\ingroup matching ///\file ///\brief Maximum matching algorithms in bipartite graphs. /// ///\note The pr_bipartite_matching.h file also contains algorithms to ///solve maximum cardinality bipartite matching problems. namespace lemon { /// \ingroup matching /// /// \brief Bipartite Max Cardinality Matching algorithm /// /// Bipartite Max Cardinality Matching algorithm. This class implements /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time /// complexity. /// /// \note In several cases the push-relabel based algorithms have /// better runtime performance than the augmenting path based ones. /// /// \see PrBipartiteMatching template class MaxBipartiteMatching { protected: typedef BpUGraph Graph; typedef typename Graph::Node Node; typedef typename Graph::ANodeIt ANodeIt; typedef typename Graph::BNodeIt BNodeIt; typedef typename Graph::UEdge UEdge; typedef typename Graph::UEdgeIt UEdgeIt; typedef typename Graph::IncEdgeIt IncEdgeIt; typedef typename BpUGraph::template ANodeMap ANodeMatchingMap; typedef typename BpUGraph::template BNodeMap BNodeMatchingMap; public: /// \brief Constructor. /// /// Constructor of the algorithm. MaxBipartiteMatching(const BpUGraph& graph) : _matching(graph), _rmatching(graph), _reached(graph), _graph(&graph) {} /// \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() or one alternative for it. /// Finally \ref start() will perform the matching computation or /// with step-by-step execution you can augment the solution. /// @{ /// \brief Initalize the data structures. /// /// It initalizes the data structures and creates an empty matching. void init() { for (ANodeIt it(*_graph); it != INVALID; ++it) { _matching.set(it, INVALID); } for (BNodeIt it(*_graph); it != INVALID; ++it) { _rmatching.set(it, INVALID); _reached.set(it, -1); } _size = 0; _phase = -1; } /// \brief Initalize the data structures. /// /// It initalizes the data structures and creates a greedy /// matching. From this matching sometimes it is faster to get /// the matching than from the initial empty matching. void greedyInit() { _size = 0; for (BNodeIt it(*_graph); it != INVALID; ++it) { _rmatching.set(it, INVALID); _reached.set(it, 0); } for (ANodeIt it(*_graph); it != INVALID; ++it) { _matching[it] = INVALID; for (IncEdgeIt jt(*_graph, it); jt != INVALID; ++jt) { if (_rmatching[_graph->bNode(jt)] == INVALID) { _matching.set(it, jt); _rmatching.set(_graph->bNode(jt), jt); _reached.set(it, -1); ++_size; break; } } } _phase = 0; } /// \brief Initalize the data structures with an initial matching. /// /// It initalizes the data structures with an initial matching. template void matchingInit(const MatchingMap& mm) { for (ANodeIt it(*_graph); it != INVALID; ++it) { _matching.set(it, INVALID); } for (BNodeIt it(*_graph); it != INVALID; ++it) { _rmatching.set(it, INVALID); _reached.set(it, 0); } _size = 0; for (UEdgeIt it(*_graph); it != INVALID; ++it) { if (mm[it]) { ++_size; _matching.set(_graph->aNode(it), it); _rmatching.set(_graph->bNode(it), it); _reached.set(it, 0); } } _phase = 0; } /// \brief Initalize the data structures with an initial matching. /// /// It initalizes the data structures with an initial matching. /// \return %True when the given map contains really a matching. template bool checkedMatchingInit(const MatchingMap& mm) { for (ANodeIt it(*_graph); it != INVALID; ++it) { _matching.set(it, INVALID); } for (BNodeIt it(*_graph); it != INVALID; ++it) { _rmatching.set(it, INVALID); _reached.set(it, 0); } _size = 0; for (UEdgeIt it(*_graph); it != INVALID; ++it) { if (mm[it]) { ++_size; if (_matching[_graph->aNode(it)] != INVALID) { return false; } _matching.set(_graph->aNode(it), it); if (_matching[_graph->bNode(it)] != INVALID) { return false; } _matching.set(_graph->bNode(it), it); _reached.set(_graph->bNode(it), -1); } } _phase = 0; return true; } private: bool _find_path(Node anode, int maxlevel, typename Graph::template BNodeMap& level) { for (IncEdgeIt it(*_graph, anode); it != INVALID; ++it) { Node bnode = _graph->bNode(it); if (level[bnode] == maxlevel) { level.set(bnode, -1); if (maxlevel == 0) { _matching.set(anode, it); _rmatching.set(bnode, it); return true; } else { Node nnode = _graph->aNode(_rmatching[bnode]); if (_find_path(nnode, maxlevel - 1, level)) { _matching.set(anode, it); _rmatching.set(bnode, it); return true; } } } } return false; } public: /// \brief An augmenting phase of the Hopcroft-Karp algorithm /// /// It runs an augmenting phase of the Hopcroft-Karp /// algorithm. This phase finds maximal edge disjoint augmenting /// paths and augments on these paths. The algorithm consists at /// most of \f$ O(\sqrt{n}) \f$ phase and one phase is \f$ O(e) /// \f$ long. bool augment() { ++_phase; typename Graph::template BNodeMap _level(*_graph, -1); typename Graph::template ANodeMap _found(*_graph, false); std::vector queue, aqueue; for (BNodeIt it(*_graph); it != INVALID; ++it) { if (_rmatching[it] == INVALID) { queue.push_back(it); _reached.set(it, _phase); _level.set(it, 0); } } bool success = false; int level = 0; while (!success && !queue.empty()) { std::vector nqueue; for (int i = 0; i < int(queue.size()); ++i) { Node bnode = queue[i]; for (IncEdgeIt jt(*_graph, bnode); jt != INVALID; ++jt) { Node anode = _graph->aNode(jt); if (_matching[anode] == INVALID) { if (!_found[anode]) { if (_find_path(anode, level, _level)) { ++_size; } _found.set(anode, true); } success = true; } else { Node nnode = _graph->bNode(_matching[anode]); if (_reached[nnode] != _phase) { _reached.set(nnode, _phase); nqueue.push_back(nnode); _level.set(nnode, level + 1); } } } } ++level; queue.swap(nqueue); } return success; } private: void _find_path_bfs(Node anode, typename Graph::template ANodeMap& pred) { while (true) { UEdge uedge = pred[anode]; Node bnode = _graph->bNode(uedge); UEdge nedge = _rmatching[bnode]; _matching.set(anode, uedge); _rmatching.set(bnode, uedge); if (nedge == INVALID) break; anode = _graph->aNode(nedge); } } public: /// \brief An augmenting phase with single path augementing /// /// This phase finds only one augmenting paths and augments on /// these paths. The algorithm consists at most of \f$ O(n) \f$ /// phase and one phase is \f$ O(e) \f$ long. bool simpleAugment() { ++_phase; typename Graph::template ANodeMap _pred(*_graph); std::vector queue, aqueue; for (BNodeIt it(*_graph); it != INVALID; ++it) { if (_rmatching[it] == INVALID) { queue.push_back(it); _reached.set(it, _phase); } } bool success = false; int level = 0; while (!success && !queue.empty()) { std::vector nqueue; for (int i = 0; i < int(queue.size()); ++i) { Node bnode = queue[i]; for (IncEdgeIt jt(*_graph, bnode); jt != INVALID; ++jt) { Node anode = _graph->aNode(jt); if (_matching[anode] == INVALID) { _pred.set(anode, jt); _find_path_bfs(anode, _pred); ++_size; return true; } else { Node nnode = _graph->bNode(_matching[anode]); if (_reached[nnode] != _phase) { _pred.set(anode, jt); _reached.set(nnode, _phase); nqueue.push_back(nnode); } } } } ++level; queue.swap(nqueue); } return success; } /// \brief Starts the algorithm. /// /// Starts the algorithm. It runs augmenting phases until the optimal /// solution reached. void start() { while (augment()) {} } /// \brief Runs the algorithm. /// /// It just initalize the algorithm and then start it. void run() { greedyInit(); start(); } /// @} /// \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 true if the given uedge is in the matching. /// /// It returns true if the given uedge is in the matching. bool matchingEdge(const UEdge& edge) const { return _matching[_graph->aNode(edge)] == edge; } /// \brief Returns the matching edge from the node. /// /// Returns the matching edge from the node. If there is not such /// edge it gives back \c INVALID. /// \note If the parameter node is a B-node then the running time is /// propotional to the degree of the node. UEdge matchingEdge(const Node& node) const { if (_graph->aNode(node)) { return _matching[node]; } else { return _rmatching[node]; } } /// \brief Set true all matching uedge in the map. /// /// Set true all matching uedge in the map. It does not change the /// value mapped to the other uedges. /// \return The number of the matching edges. template int quickMatching(MatchingMap& mm) const { for (ANodeIt it(*_graph); it != INVALID; ++it) { if (_matching[it] != INVALID) { mm.set(_matching[it], true); } } return _size; } /// \brief Set true all matching uedge in the map and the others to false. /// /// Set true all matching uedge in the map and the others to false. /// \return The number of the matching edges. template int matching(MatchingMap& mm) const { for (UEdgeIt it(*_graph); it != INVALID; ++it) { mm.set(it, it == _matching[_graph->aNode(it)]); } return _size; } ///Gives back the matching in an ANodeMap. ///Gives back the matching in an ANodeMap. The parameter should ///be a write ANodeMap of UEdge values. ///\return The number of the matching edges. template int aMatching(MatchingMap& mm) const { for (ANodeIt it(*_graph); it != INVALID; ++it) { mm.set(it, _matching[it]); } return _size; } ///Gives back the matching in a BNodeMap. ///Gives back the matching in a BNodeMap. The parameter should ///be a write BNodeMap of UEdge values. ///\return The number of the matching edges. template int bMatching(MatchingMap& mm) const { for (BNodeIt it(*_graph); it != INVALID; ++it) { mm.set(it, _rmatching[it]); } return _size; } /// \brief Returns a minimum covering of the nodes. /// /// The minimum covering set problem is the dual solution of the /// maximum bipartite matching. It provides a solution for this /// problem what is proof of the optimality of the matching. /// \return The size of the cover set. template int coverSet(CoverMap& covering) const { int size = 0; for (ANodeIt it(*_graph); it != INVALID; ++it) { bool cn = _matching[it] != INVALID && _reached[_graph->bNode(_matching[it])] == _phase; covering.set(it, cn); if (cn) ++size; } for (BNodeIt it(*_graph); it != INVALID; ++it) { bool cn = _reached[it] != _phase; covering.set(it, cn); if (cn) ++size; } return size; } /// \brief Gives back a barrier on the A-nodes /// /// The barrier is s subset of the nodes on the same side of the /// graph, which size minus its neighbours is exactly the /// unmatched nodes on the A-side. /// \retval barrier A WriteMap on the ANodes with bool value. template void aBarrier(BarrierMap& barrier) const { for (ANodeIt it(*_graph); it != INVALID; ++it) { barrier.set(it, _matching[it] == INVALID || _reached[_graph->bNode(_matching[it])] != _phase); } } /// \brief Gives back a barrier on the B-nodes /// /// The barrier is s subset of the nodes on the same side of the /// graph, which size minus its neighbours is exactly the /// unmatched nodes on the B-side. /// \retval barrier A WriteMap on the BNodes with bool value. template void bBarrier(BarrierMap& barrier) const { for (BNodeIt it(*_graph); it != INVALID; ++it) { barrier.set(it, _reached[it] == _phase); } } /// \brief Gives back the number of the matching edges. /// /// Gives back the number of the matching edges. int matchingSize() const { return _size; } /// @} private: typename BpUGraph::template ANodeMap _matching; typename BpUGraph::template BNodeMap _rmatching; typename BpUGraph::template BNodeMap _reached; int _phase; const Graph *_graph; int _size; }; /// \ingroup matching /// /// \brief Maximum cardinality bipartite matching /// /// This function calculates the maximum cardinality matching /// in a bipartite graph. It gives back the matching in an undirected /// edge map. /// /// \param graph The bipartite graph. /// \return The size of the matching. template int maxBipartiteMatching(const BpUGraph& graph) { MaxBipartiteMatching bpmatching(graph); bpmatching.run(); return bpmatching.matchingSize(); } /// \ingroup matching /// /// \brief Maximum cardinality bipartite matching /// /// This function calculates the maximum cardinality matching /// in a bipartite graph. It gives back the matching in an undirected /// edge map. /// /// \param graph The bipartite graph. /// \retval matching The ANodeMap of UEdges which will be set to covered /// matching undirected edge. /// \return The size of the matching. template int maxBipartiteMatching(const BpUGraph& graph, MatchingMap& matching) { MaxBipartiteMatching bpmatching(graph); bpmatching.run(); bpmatching.aMatching(matching); return bpmatching.matchingSize(); } /// \ingroup matching /// /// \brief Maximum cardinality bipartite matching /// /// This function calculates the maximum cardinality matching /// in a bipartite graph. It gives back the matching in an undirected /// edge map. /// /// \param graph The bipartite graph. /// \retval matching The ANodeMap of UEdges which will be set to covered /// matching undirected edge. /// \retval barrier The BNodeMap of bools which will be set to a barrier /// of the BNode-set. /// \return The size of the matching. template int maxBipartiteMatching(const BpUGraph& graph, MatchingMap& matching, BarrierMap& barrier) { MaxBipartiteMatching bpmatching(graph); bpmatching.run(); bpmatching.aMatching(matching); bpmatching.bBarrier(barrier); return bpmatching.matchingSize(); } /// \brief Default traits class for weighted bipartite matching algoritms. /// /// Default traits class for weighted bipartite matching algoritms. /// \param _BpUGraph The bipartite undirected graph type. /// \param _WeightMap Type of weight map. template struct WeightedBipartiteMatchingDefaultTraits { /// \brief The type of the weight of the undirected edges. typedef typename _WeightMap::Value Value; /// The undirected bipartite graph type the algorithm runs on. typedef _BpUGraph BpUGraph; /// The map of the edges weights typedef _WeightMap WeightMap; /// \brief The cross reference type used by heap. /// /// The cross reference type used by heap. /// Usually it is \c Graph::NodeMap. typedef typename BpUGraph::template NodeMap HeapCrossRef; /// \brief Instantiates a HeapCrossRef. /// /// This function instantiates a \ref HeapCrossRef. /// \param graph is the graph, to which we would like to define the /// HeapCrossRef. static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) { return new HeapCrossRef(graph); } /// \brief The heap type used by weighted matching algorithms. /// /// The heap type used by weighted matching algorithms. It should /// minimize the priorities and the heap's key type is the graph's /// anode graph's node. /// /// \sa BinHeap typedef BinHeap Heap; /// \brief Instantiates a Heap. /// /// This function instantiates a \ref Heap. /// \param crossref The cross reference of the heap. static Heap *createHeap(HeapCrossRef& crossref) { return new Heap(crossref); } }; /// \ingroup matching /// /// \brief Bipartite Max Weighted Matching algorithm /// /// This class implements the bipartite Max Weighted Matching /// algorithm. It uses the successive shortest path algorithm to /// calculate the maximum weighted matching in the bipartite /// graph. The algorithm can be used also to calculate the maximum /// cardinality maximum weighted matching. The time complexity /// of the algorithm is \f$ O(ne\log(n)) \f$ with the default binary /// heap implementation but this can be improved to /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps. /// /// The algorithm also provides a potential function on the nodes /// which a dual solution of the matching algorithm and it can be /// used to proof the optimality of the given pimal solution. #ifdef DOXYGEN template #else template , typename _Traits = WeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> > #endif class MaxWeightedBipartiteMatching { public: typedef _Traits Traits; typedef typename Traits::BpUGraph BpUGraph; typedef typename Traits::WeightMap WeightMap; typedef typename Traits::Value Value; protected: typedef typename Traits::HeapCrossRef HeapCrossRef; typedef typename Traits::Heap Heap; typedef typename BpUGraph::Node Node; typedef typename BpUGraph::ANodeIt ANodeIt; typedef typename BpUGraph::BNodeIt BNodeIt; typedef typename BpUGraph::UEdge UEdge; typedef typename BpUGraph::UEdgeIt UEdgeIt; typedef typename BpUGraph::IncEdgeIt IncEdgeIt; typedef typename BpUGraph::template ANodeMap ANodeMatchingMap; typedef typename BpUGraph::template BNodeMap BNodeMatchingMap; typedef typename BpUGraph::template ANodeMap ANodePotentialMap; typedef typename BpUGraph::template BNodeMap BNodePotentialMap; public: /// \brief \ref Exception for uninitialized parameters. /// /// This error represents problems in the initialization /// of the parameters of the algorithms. class UninitializedParameter : public lemon::UninitializedParameter { public: virtual const char* what() const throw() { return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter"; } }; ///\name Named template parameters ///@{ template struct DefHeapTraits : public Traits { typedef CR HeapCrossRef; typedef H Heap; static HeapCrossRef *createHeapCrossRef(const BpUGraph &) { throw UninitializedParameter(); } static Heap *createHeap(HeapCrossRef &) { throw UninitializedParameter(); } }; /// \brief \ref named-templ-param "Named parameter" for setting heap /// and cross reference type /// /// \ref named-templ-param "Named parameter" for setting heap and cross /// reference type template > struct DefHeap : public MaxWeightedBipartiteMatching > { typedef MaxWeightedBipartiteMatching > Create; }; template struct DefStandardHeapTraits : public Traits { typedef CR HeapCrossRef; typedef H Heap; static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) { return new HeapCrossRef(graph); } static Heap *createHeap(HeapCrossRef &crossref) { return new Heap(crossref); } }; /// \brief \ref named-templ-param "Named parameter" for setting heap and /// cross reference type with automatic allocation /// /// \ref named-templ-param "Named parameter" for setting heap and cross /// reference type. It can allocate the heap and the cross reference /// object if the cross reference's constructor waits for the graph as /// parameter and the heap's constructor waits for the cross reference. template > struct DefStandardHeap : public MaxWeightedBipartiteMatching > { typedef MaxWeightedBipartiteMatching > Create; }; ///@} /// \brief Constructor. /// /// Constructor of the algorithm. MaxWeightedBipartiteMatching(const BpUGraph& _graph, const WeightMap& _weight) : graph(&_graph), weight(&_weight), anode_matching(_graph), bnode_matching(_graph), anode_potential(_graph), bnode_potential(_graph), _heap_cross_ref(0), local_heap_cross_ref(false), _heap(0), local_heap(0) {} /// \brief Destructor. /// /// Destructor of the algorithm. ~MaxWeightedBipartiteMatching() { destroyStructures(); } /// \brief Sets the heap and the cross reference used by algorithm. /// /// Sets the heap and the cross reference used by algorithm. /// If you don't use this function before calling \ref run(), /// it will allocate one. The destuctor deallocates this /// automatically allocated map, of course. /// \return \c (*this) MaxWeightedBipartiteMatching& heap(Heap& hp, HeapCrossRef &cr) { if(local_heap_cross_ref) { delete _heap_cross_ref; local_heap_cross_ref = false; } _heap_cross_ref = &cr; if(local_heap) { delete _heap; local_heap = false; } _heap = &hp; return *this; } /// \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() or one alternative for it. /// Finally \ref start() will perform the matching computation or /// with step-by-step execution you can augment the solution. /// @{ /// \brief Initalize the data structures. /// /// It initalizes the data structures and creates an empty matching. void init() { initStructures(); for (ANodeIt it(*graph); it != INVALID; ++it) { anode_matching[it] = INVALID; anode_potential[it] = 0; } for (BNodeIt it(*graph); it != INVALID; ++it) { bnode_matching[it] = INVALID; bnode_potential[it] = 0; for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) { if ((*weight)[jt] > bnode_potential[it]) { bnode_potential[it] = (*weight)[jt]; } } } matching_value = 0; matching_size = 0; } /// \brief An augmenting phase of the weighted matching algorithm /// /// It runs an augmenting phase of the weighted matching /// algorithm. This phase finds the best augmenting path and /// augments only on this paths. /// /// The algorithm consists at most /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long /// with binary heap. /// \param decrease If the given parameter true the matching value /// can be decreased in the augmenting phase. If we would like /// to calculate the maximum cardinality maximum weighted matching /// then we should let the algorithm to decrease the matching /// value in order to increase the number of the matching edges. bool augment(bool decrease = false) { typename BpUGraph::template BNodeMap bdist(*graph); typename BpUGraph::template BNodeMap bpred(*graph, INVALID); Node bestNode = INVALID; Value bestValue = 0; _heap->clear(); for (ANodeIt it(*graph); it != INVALID; ++it) { (*_heap_cross_ref)[it] = Heap::PRE_HEAP; } for (ANodeIt it(*graph); it != INVALID; ++it) { if (anode_matching[it] == INVALID) { _heap->push(it, 0); } } Value bdistMax = 0; while (!_heap->empty()) { Node anode = _heap->top(); Value avalue = _heap->prio(); _heap->pop(); for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) { if (jt == anode_matching[anode]) continue; Node bnode = graph->bNode(jt); Value bvalue = avalue - (*weight)[jt] + anode_potential[anode] + bnode_potential[bnode]; if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) { bdist[bnode] = bvalue; bpred[bnode] = jt; } if (bvalue > bdistMax) { bdistMax = bvalue; } if (bnode_matching[bnode] != INVALID) { Node newanode = graph->aNode(bnode_matching[bnode]); switch (_heap->state(newanode)) { case Heap::PRE_HEAP: _heap->push(newanode, bvalue); break; case Heap::IN_HEAP: if (bvalue < (*_heap)[newanode]) { _heap->decrease(newanode, bvalue); } break; case Heap::POST_HEAP: break; } } else { if (bestNode == INVALID || bnode_potential[bnode] - bvalue > bestValue) { bestValue = bnode_potential[bnode] - bvalue; bestNode = bnode; } } } } if (bestNode == INVALID || (!decrease && bestValue < 0)) { return false; } matching_value += bestValue; ++matching_size; for (BNodeIt it(*graph); it != INVALID; ++it) { if (bpred[it] != INVALID) { bnode_potential[it] -= bdist[it]; } else { bnode_potential[it] -= bdistMax; } } for (ANodeIt it(*graph); it != INVALID; ++it) { if (anode_matching[it] != INVALID) { Node bnode = graph->bNode(anode_matching[it]); if (bpred[bnode] != INVALID) { anode_potential[it] += bdist[bnode]; } else { anode_potential[it] += bdistMax; } } } while (bestNode != INVALID) { UEdge uedge = bpred[bestNode]; Node anode = graph->aNode(uedge); bnode_matching[bestNode] = uedge; if (anode_matching[anode] != INVALID) { bestNode = graph->bNode(anode_matching[anode]); } else { bestNode = INVALID; } anode_matching[anode] = uedge; } return true; } /// \brief Starts the algorithm. /// /// Starts the algorithm. It runs augmenting phases until the /// optimal solution reached. /// /// \param maxCardinality If the given value is true it will /// calculate the maximum cardinality maximum matching instead of /// the maximum matching. void start(bool maxCardinality = false) { while (augment(maxCardinality)) {} } /// \brief Runs the algorithm. /// /// It just initalize the algorithm and then start it. /// /// \param maxCardinality If the given value is true it will /// calculate the maximum cardinality maximum matching instead of /// the maximum matching. void run(bool maxCardinality = false) { init(); start(maxCardinality); } /// @} /// \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 Gives back the potential in the NodeMap /// /// Gives back the potential in the NodeMap. The matching is optimal /// with the current number of edges if \f$ \pi(a) + \pi(b) - w(ab) = 0 \f$ /// for each matching edges and \f$ \pi(a) + \pi(b) - w(ab) \ge 0 \f$ /// for each edges. template void potential(PotentialMap& pt) const { for (ANodeIt it(*graph); it != INVALID; ++it) { pt.set(it, anode_potential[it]); } for (BNodeIt it(*graph); it != INVALID; ++it) { pt.set(it, bnode_potential[it]); } } /// \brief Set true all matching uedge in the map. /// /// Set true all matching uedge in the map. It does not change the /// value mapped to the other uedges. /// \return The number of the matching edges. template int quickMatching(MatchingMap& mm) const { for (ANodeIt it(*graph); it != INVALID; ++it) { if (anode_matching[it] != INVALID) { mm.set(anode_matching[it], true); } } return matching_size; } /// \brief Set true all matching uedge in the map and the others to false. /// /// Set true all matching uedge in the map and the others to false. /// \return The number of the matching edges. template int matching(MatchingMap& mm) const { for (UEdgeIt it(*graph); it != INVALID; ++it) { mm.set(it, it == anode_matching[graph->aNode(it)]); } return matching_size; } ///Gives back the matching in an ANodeMap. ///Gives back the matching in an ANodeMap. The parameter should ///be a write ANodeMap of UEdge values. ///\return The number of the matching edges. template int aMatching(MatchingMap& mm) const { for (ANodeIt it(*graph); it != INVALID; ++it) { mm.set(it, anode_matching[it]); } return matching_size; } ///Gives back the matching in a BNodeMap. ///Gives back the matching in a BNodeMap. The parameter should ///be a write BNodeMap of UEdge values. ///\return The number of the matching edges. template int bMatching(MatchingMap& mm) const { for (BNodeIt it(*graph); it != INVALID; ++it) { mm.set(it, bnode_matching[it]); } return matching_size; } /// \brief Return true if the given uedge is in the matching. /// /// It returns true if the given uedge is in the matching. bool matchingEdge(const UEdge& edge) const { return anode_matching[graph->aNode(edge)] == edge; } /// \brief Returns the matching edge from the node. /// /// Returns the matching edge from the node. If there is not such /// edge it gives back \c INVALID. UEdge matchingEdge(const Node& node) const { if (graph->aNode(node)) { return anode_matching[node]; } else { return bnode_matching[node]; } } /// \brief Gives back the sum of weights of the matching edges. /// /// Gives back the sum of weights of the matching edges. Value matchingValue() const { return matching_value; } /// \brief Gives back the number of the matching edges. /// /// Gives back the number of the matching edges. int matchingSize() const { return matching_size; } /// @} private: void initStructures() { if (!_heap_cross_ref) { local_heap_cross_ref = true; _heap_cross_ref = Traits::createHeapCrossRef(*graph); } if (!_heap) { local_heap = true; _heap = Traits::createHeap(*_heap_cross_ref); } } void destroyStructures() { if (local_heap_cross_ref) delete _heap_cross_ref; if (local_heap) delete _heap; } private: const BpUGraph *graph; const WeightMap* weight; ANodeMatchingMap anode_matching; BNodeMatchingMap bnode_matching; ANodePotentialMap anode_potential; BNodePotentialMap bnode_potential; Value matching_value; int matching_size; HeapCrossRef *_heap_cross_ref; bool local_heap_cross_ref; Heap *_heap; bool local_heap; }; /// \ingroup matching /// /// \brief Maximum weighted bipartite matching /// /// This function calculates the maximum weighted matching /// in a bipartite graph. It gives back the matching in an undirected /// edge map. /// /// \param graph The bipartite graph. /// \param weight The undirected edge map which contains the weights. /// \retval matching The undirected edge map which will be set to /// the matching. /// \return The value of the matching. template typename WeightMap::Value maxWeightedBipartiteMatching(const BpUGraph& graph, const WeightMap& weight, MatchingMap& matching) { MaxWeightedBipartiteMatching bpmatching(graph, weight); bpmatching.run(); bpmatching.matching(matching); return bpmatching.matchingValue(); } /// \ingroup matching /// /// \brief Maximum weighted maximum cardinality bipartite matching /// /// This function calculates the maximum weighted of the maximum cardinality /// matchings of a bipartite graph. It gives back the matching in an /// undirected edge map. /// /// \param graph The bipartite graph. /// \param weight The undirected edge map which contains the weights. /// \retval matching The undirected edge map which will be set to /// the matching. /// \return The value of the matching. template typename WeightMap::Value maxWeightedMaxBipartiteMatching(const BpUGraph& graph, const WeightMap& weight, MatchingMap& matching) { MaxWeightedBipartiteMatching bpmatching(graph, weight); bpmatching.run(true); bpmatching.matching(matching); return bpmatching.matchingValue(); } /// \brief Default traits class for minimum cost bipartite matching /// algoritms. /// /// Default traits class for minimum cost bipartite matching /// algoritms. /// /// \param _BpUGraph The bipartite undirected graph /// type. /// /// \param _CostMap Type of cost map. template struct MinCostMaxBipartiteMatchingDefaultTraits { /// \brief The type of the cost of the undirected edges. typedef typename _CostMap::Value Value; /// The undirected bipartite graph type the algorithm runs on. typedef _BpUGraph BpUGraph; /// The map of the edges costs typedef _CostMap CostMap; /// \brief The cross reference type used by heap. /// /// The cross reference type used by heap. /// Usually it is \c Graph::NodeMap. typedef typename BpUGraph::template NodeMap HeapCrossRef; /// \brief Instantiates a HeapCrossRef. /// /// This function instantiates a \ref HeapCrossRef. /// \param graph is the graph, to which we would like to define the /// HeapCrossRef. static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) { return new HeapCrossRef(graph); } /// \brief The heap type used by costed matching algorithms. /// /// The heap type used by costed matching algorithms. It should /// minimize the priorities and the heap's key type is the graph's /// anode graph's node. /// /// \sa BinHeap typedef BinHeap Heap; /// \brief Instantiates a Heap. /// /// This function instantiates a \ref Heap. /// \param crossref The cross reference of the heap. static Heap *createHeap(HeapCrossRef& crossref) { return new Heap(crossref); } }; /// \ingroup matching /// /// \brief Bipartite Min Cost Matching algorithm /// /// This class implements the bipartite Min Cost Matching algorithm. /// It uses the successive shortest path algorithm to calculate the /// minimum cost maximum matching in the bipartite graph. The time /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the /// default binary heap implementation but this can be improved to /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps. /// /// The algorithm also provides a potential function on the nodes /// which a dual solution of the matching algorithm and it can be /// used to proof the optimality of the given pimal solution. #ifdef DOXYGEN template #else template , typename _Traits = MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> > #endif class MinCostMaxBipartiteMatching { public: typedef _Traits Traits; typedef typename Traits::BpUGraph BpUGraph; typedef typename Traits::CostMap CostMap; typedef typename Traits::Value Value; protected: typedef typename Traits::HeapCrossRef HeapCrossRef; typedef typename Traits::Heap Heap; typedef typename BpUGraph::Node Node; typedef typename BpUGraph::ANodeIt ANodeIt; typedef typename BpUGraph::BNodeIt BNodeIt; typedef typename BpUGraph::UEdge UEdge; typedef typename BpUGraph::UEdgeIt UEdgeIt; typedef typename BpUGraph::IncEdgeIt IncEdgeIt; typedef typename BpUGraph::template ANodeMap ANodeMatchingMap; typedef typename BpUGraph::template BNodeMap BNodeMatchingMap; typedef typename BpUGraph::template ANodeMap ANodePotentialMap; typedef typename BpUGraph::template BNodeMap BNodePotentialMap; public: /// \brief \ref Exception for uninitialized parameters. /// /// This error represents problems in the initialization /// of the parameters of the algorithms. class UninitializedParameter : public lemon::UninitializedParameter { public: virtual const char* what() const throw() { return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter"; } }; ///\name Named template parameters ///@{ template struct DefHeapTraits : public Traits { typedef CR HeapCrossRef; typedef H Heap; static HeapCrossRef *createHeapCrossRef(const BpUGraph &) { throw UninitializedParameter(); } static Heap *createHeap(HeapCrossRef &) { throw UninitializedParameter(); } }; /// \brief \ref named-templ-param "Named parameter" for setting heap /// and cross reference type /// /// \ref named-templ-param "Named parameter" for setting heap and cross /// reference type template > struct DefHeap : public MinCostMaxBipartiteMatching > { typedef MinCostMaxBipartiteMatching > Create; }; template struct DefStandardHeapTraits : public Traits { typedef CR HeapCrossRef; typedef H Heap; static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) { return new HeapCrossRef(graph); } static Heap *createHeap(HeapCrossRef &crossref) { return new Heap(crossref); } }; /// \brief \ref named-templ-param "Named parameter" for setting heap and /// cross reference type with automatic allocation /// /// \ref named-templ-param "Named parameter" for setting heap and cross /// reference type. It can allocate the heap and the cross reference /// object if the cross reference's constructor waits for the graph as /// parameter and the heap's constructor waits for the cross reference. template > struct DefStandardHeap : public MinCostMaxBipartiteMatching > { typedef MinCostMaxBipartiteMatching > Create; }; ///@} /// \brief Constructor. /// /// Constructor of the algorithm. MinCostMaxBipartiteMatching(const BpUGraph& _graph, const CostMap& _cost) : graph(&_graph), cost(&_cost), anode_matching(_graph), bnode_matching(_graph), anode_potential(_graph), bnode_potential(_graph), _heap_cross_ref(0), local_heap_cross_ref(false), _heap(0), local_heap(0) {} /// \brief Destructor. /// /// Destructor of the algorithm. ~MinCostMaxBipartiteMatching() { destroyStructures(); } /// \brief Sets the heap and the cross reference used by algorithm. /// /// Sets the heap and the cross reference used by algorithm. /// If you don't use this function before calling \ref run(), /// it will allocate one. The destuctor deallocates this /// automatically allocated map, of course. /// \return \c (*this) MinCostMaxBipartiteMatching& heap(Heap& hp, HeapCrossRef &cr) { if(local_heap_cross_ref) { delete _heap_cross_ref; local_heap_cross_ref = false; } _heap_cross_ref = &cr; if(local_heap) { delete _heap; local_heap = false; } _heap = &hp; return *this; } /// \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() or one alternative for it. /// Finally \ref start() will perform the matching computation or /// with step-by-step execution you can augment the solution. /// @{ /// \brief Initalize the data structures. /// /// It initalizes the data structures and creates an empty matching. void init() { initStructures(); for (ANodeIt it(*graph); it != INVALID; ++it) { anode_matching[it] = INVALID; anode_potential[it] = 0; } for (BNodeIt it(*graph); it != INVALID; ++it) { bnode_matching[it] = INVALID; bnode_potential[it] = 0; } matching_cost = 0; matching_size = 0; } /// \brief An augmenting phase of the costed matching algorithm /// /// It runs an augmenting phase of the matching algorithm. The /// phase finds the best augmenting path and augments only on this /// paths. /// /// The algorithm consists at most /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long /// with binary heap. bool augment() { typename BpUGraph::template BNodeMap bdist(*graph); typename BpUGraph::template BNodeMap bpred(*graph, INVALID); Node bestNode = INVALID; Value bestValue = 0; _heap->clear(); for (ANodeIt it(*graph); it != INVALID; ++it) { (*_heap_cross_ref)[it] = Heap::PRE_HEAP; } for (ANodeIt it(*graph); it != INVALID; ++it) { if (anode_matching[it] == INVALID) { _heap->push(it, 0); } } Value bdistMax = 0; while (!_heap->empty()) { Node anode = _heap->top(); Value avalue = _heap->prio(); _heap->pop(); for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) { if (jt == anode_matching[anode]) continue; Node bnode = graph->bNode(jt); Value bvalue = avalue + (*cost)[jt] + anode_potential[anode] - bnode_potential[bnode]; if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) { bdist[bnode] = bvalue; bpred[bnode] = jt; } if (bvalue > bdistMax) { bdistMax = bvalue; } if (bnode_matching[bnode] != INVALID) { Node newanode = graph->aNode(bnode_matching[bnode]); switch (_heap->state(newanode)) { case Heap::PRE_HEAP: _heap->push(newanode, bvalue); break; case Heap::IN_HEAP: if (bvalue < (*_heap)[newanode]) { _heap->decrease(newanode, bvalue); } break; case Heap::POST_HEAP: break; } } else { if (bestNode == INVALID || bvalue + bnode_potential[bnode] < bestValue) { bestValue = bvalue + bnode_potential[bnode]; bestNode = bnode; } } } } if (bestNode == INVALID) { return false; } matching_cost += bestValue; ++matching_size; for (BNodeIt it(*graph); it != INVALID; ++it) { if (bpred[it] != INVALID) { bnode_potential[it] += bdist[it]; } else { bnode_potential[it] += bdistMax; } } for (ANodeIt it(*graph); it != INVALID; ++it) { if (anode_matching[it] != INVALID) { Node bnode = graph->bNode(anode_matching[it]); if (bpred[bnode] != INVALID) { anode_potential[it] += bdist[bnode]; } else { anode_potential[it] += bdistMax; } } } while (bestNode != INVALID) { UEdge uedge = bpred[bestNode]; Node anode = graph->aNode(uedge); bnode_matching[bestNode] = uedge; if (anode_matching[anode] != INVALID) { bestNode = graph->bNode(anode_matching[anode]); } else { bestNode = INVALID; } anode_matching[anode] = uedge; } return true; } /// \brief Starts the algorithm. /// /// Starts the algorithm. It runs augmenting phases until the /// optimal solution reached. void start() { while (augment()) {} } /// \brief Runs the algorithm. /// /// It just initalize the algorithm and then start it. void run() { init(); start(); } /// @} /// \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 Gives back the potential in the NodeMap /// /// Gives back the potential in the NodeMap. The matching is optimal /// with the current number of edges if \f$ \pi(a) + \pi(b) - w(ab) = 0 \f$ /// for each matching edges and \f$ \pi(a) + \pi(b) - w(ab) \ge 0 \f$ /// for each edges. template void potential(PotentialMap& pt) const { for (ANodeIt it(*graph); it != INVALID; ++it) { pt.set(it, anode_potential[it]); } for (BNodeIt it(*graph); it != INVALID; ++it) { pt.set(it, bnode_potential[it]); } } /// \brief Set true all matching uedge in the map. /// /// Set true all matching uedge in the map. It does not change the /// value mapped to the other uedges. /// \return The number of the matching edges. template int quickMatching(MatchingMap& mm) const { for (ANodeIt it(*graph); it != INVALID; ++it) { if (anode_matching[it] != INVALID) { mm.set(anode_matching[it], true); } } return matching_size; } /// \brief Set true all matching uedge in the map and the others to false. /// /// Set true all matching uedge in the map and the others to false. /// \return The number of the matching edges. template int matching(MatchingMap& mm) const { for (UEdgeIt it(*graph); it != INVALID; ++it) { mm.set(it, it == anode_matching[graph->aNode(it)]); } return matching_size; } /// \brief Gives back the matching in an ANodeMap. /// /// Gives back the matching in an ANodeMap. The parameter should /// be a write ANodeMap of UEdge values. /// \return The number of the matching edges. template int aMatching(MatchingMap& mm) const { for (ANodeIt it(*graph); it != INVALID; ++it) { mm.set(it, anode_matching[it]); } return matching_size; } /// \brief Gives back the matching in a BNodeMap. /// /// Gives back the matching in a BNodeMap. The parameter should /// be a write BNodeMap of UEdge values. /// \return The number of the matching edges. template int bMatching(MatchingMap& mm) const { for (BNodeIt it(*graph); it != INVALID; ++it) { mm.set(it, bnode_matching[it]); } return matching_size; } /// \brief Return true if the given uedge is in the matching. /// /// It returns true if the given uedge is in the matching. bool matchingEdge(const UEdge& edge) const { return anode_matching[graph->aNode(edge)] == edge; } /// \brief Returns the matching edge from the node. /// /// Returns the matching edge from the node. If there is not such /// edge it gives back \c INVALID. UEdge matchingEdge(const Node& node) const { if (graph->aNode(node)) { return anode_matching[node]; } else { return bnode_matching[node]; } } /// \brief Gives back the sum of costs of the matching edges. /// /// Gives back the sum of costs of the matching edges. Value matchingCost() const { return matching_cost; } /// \brief Gives back the number of the matching edges. /// /// Gives back the number of the matching edges. int matchingSize() const { return matching_size; } /// @} private: void initStructures() { if (!_heap_cross_ref) { local_heap_cross_ref = true; _heap_cross_ref = Traits::createHeapCrossRef(*graph); } if (!_heap) { local_heap = true; _heap = Traits::createHeap(*_heap_cross_ref); } } void destroyStructures() { if (local_heap_cross_ref) delete _heap_cross_ref; if (local_heap) delete _heap; } private: const BpUGraph *graph; const CostMap* cost; ANodeMatchingMap anode_matching; BNodeMatchingMap bnode_matching; ANodePotentialMap anode_potential; BNodePotentialMap bnode_potential; Value matching_cost; int matching_size; HeapCrossRef *_heap_cross_ref; bool local_heap_cross_ref; Heap *_heap; bool local_heap; }; /// \ingroup matching /// /// \brief Minimum cost maximum cardinality bipartite matching /// /// This function calculates the maximum cardinality matching with /// minimum cost of a bipartite graph. It gives back the matching in /// an undirected edge map. /// /// \param graph The bipartite graph. /// \param cost The undirected edge map which contains the costs. /// \retval matching The undirected edge map which will be set to /// the matching. /// \return The cost of the matching. template typename CostMap::Value minCostMaxBipartiteMatching(const BpUGraph& graph, const CostMap& cost, MatchingMap& matching) { MinCostMaxBipartiteMatching bpmatching(graph, cost); bpmatching.run(); bpmatching.matching(matching); return bpmatching.matchingCost(); } } #endif