deba@2040: /* -*- C++ -*- deba@2040: * deba@2040: * This file is a part of LEMON, a generic C++ optimization library deba@2040: * deba@2040: * Copyright (C) 2003-2006 deba@2040: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport deba@2040: * (Egervary Research Group on Combinatorial Optimization, EGRES). deba@2040: * deba@2040: * Permission to use, modify and distribute this software is granted deba@2040: * provided that this copyright notice appears in all copies. For deba@2040: * precise terms see the accompanying LICENSE file. deba@2040: * deba@2040: * This software is provided "AS IS" with no warranty of any kind, deba@2040: * express or implied, and with no claim as to its suitability for any deba@2040: * purpose. deba@2040: * deba@2040: */ deba@2040: deba@2040: #ifndef LEMON_BIPARTITE_MATCHING deba@2040: #define LEMON_BIPARTITE_MATCHING deba@2040: deba@2051: #include deba@2051: deba@2051: #include deba@2051: #include deba@2040: deba@2040: #include deba@2040: deba@2040: ///\ingroup matching deba@2040: ///\file deba@2040: ///\brief Maximum matching algorithms in bipartite graphs. deba@2040: deba@2040: namespace lemon { deba@2040: deba@2040: /// \ingroup matching deba@2040: /// deba@2040: /// \brief Bipartite Max Cardinality Matching algorithm deba@2040: /// deba@2040: /// Bipartite Max Cardinality Matching algorithm. This class implements deba@2051: /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time deba@2040: /// complexity. deba@2040: template deba@2040: class MaxBipartiteMatching { deba@2040: protected: deba@2040: deba@2040: typedef BpUGraph Graph; deba@2040: deba@2040: typedef typename Graph::Node Node; deba@2040: typedef typename Graph::ANodeIt ANodeIt; deba@2040: typedef typename Graph::BNodeIt BNodeIt; deba@2040: typedef typename Graph::UEdge UEdge; deba@2040: typedef typename Graph::UEdgeIt UEdgeIt; deba@2040: typedef typename Graph::IncEdgeIt IncEdgeIt; deba@2040: deba@2040: typedef typename BpUGraph::template ANodeMap ANodeMatchingMap; deba@2040: typedef typename BpUGraph::template BNodeMap BNodeMatchingMap; deba@2040: deba@2040: deba@2040: public: deba@2040: deba@2040: /// \brief Constructor. deba@2040: /// deba@2040: /// Constructor of the algorithm. deba@2040: MaxBipartiteMatching(const BpUGraph& _graph) deba@2040: : anode_matching(_graph), bnode_matching(_graph), graph(&_graph) {} deba@2040: deba@2040: /// \name Execution control deba@2040: /// The simplest way to execute the algorithm is to use deba@2040: /// one of the member functions called \c run(). deba@2040: /// \n deba@2040: /// If you need more control on the execution, deba@2040: /// first you must call \ref init() or one alternative for it. deba@2040: /// Finally \ref start() will perform the matching computation or deba@2040: /// with step-by-step execution you can augment the solution. deba@2040: deba@2040: /// @{ deba@2040: deba@2040: /// \brief Initalize the data structures. deba@2040: /// deba@2040: /// It initalizes the data structures and creates an empty matching. deba@2040: void init() { deba@2040: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2040: anode_matching[it] = INVALID; deba@2040: } deba@2040: for (BNodeIt it(*graph); it != INVALID; ++it) { deba@2040: bnode_matching[it] = INVALID; deba@2040: } deba@2051: matching_size = 0; deba@2040: } deba@2040: deba@2040: /// \brief Initalize the data structures. deba@2040: /// deba@2040: /// It initalizes the data structures and creates a greedy deba@2040: /// matching. From this matching sometimes it is faster to get deba@2040: /// the matching than from the initial empty matching. deba@2040: void greedyInit() { deba@2051: matching_size = 0; deba@2040: for (BNodeIt it(*graph); it != INVALID; ++it) { deba@2040: bnode_matching[it] = INVALID; deba@2040: } deba@2040: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2040: anode_matching[it] = INVALID; deba@2040: for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) { deba@2040: if (bnode_matching[graph->bNode(jt)] == INVALID) { deba@2040: anode_matching[it] = jt; deba@2040: bnode_matching[graph->bNode(jt)] = jt; deba@2051: ++matching_size; deba@2040: break; deba@2040: } deba@2040: } deba@2040: } deba@2040: } deba@2040: deba@2040: /// \brief Initalize the data structures with an initial matching. deba@2040: /// deba@2040: /// It initalizes the data structures with an initial matching. deba@2040: template deba@2040: void matchingInit(const MatchingMap& matching) { deba@2040: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2040: anode_matching[it] = INVALID; deba@2040: } deba@2040: for (BNodeIt it(*graph); it != INVALID; ++it) { deba@2040: bnode_matching[it] = INVALID; deba@2040: } deba@2051: matching_size = 0; deba@2040: for (UEdgeIt it(*graph); it != INVALID; ++it) { deba@2040: if (matching[it]) { deba@2051: ++matching_size; deba@2040: anode_matching[graph->aNode(it)] = it; deba@2040: bnode_matching[graph->bNode(it)] = it; deba@2040: } deba@2040: } deba@2040: } deba@2040: deba@2040: /// \brief Initalize the data structures with an initial matching. deba@2040: /// deba@2040: /// It initalizes the data structures with an initial matching. deba@2040: /// \return %True when the given map contains really a matching. deba@2040: template deba@2040: void checkedMatchingInit(const MatchingMap& matching) { deba@2040: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2040: anode_matching[it] = INVALID; deba@2040: } deba@2040: for (BNodeIt it(*graph); it != INVALID; ++it) { deba@2040: bnode_matching[it] = INVALID; deba@2040: } deba@2051: matching_size = 0; deba@2040: for (UEdgeIt it(*graph); it != INVALID; ++it) { deba@2040: if (matching[it]) { deba@2051: ++matching_size; deba@2040: if (anode_matching[graph->aNode(it)] != INVALID) { deba@2040: return false; deba@2040: } deba@2040: anode_matching[graph->aNode(it)] = it; deba@2040: if (bnode_matching[graph->aNode(it)] != INVALID) { deba@2040: return false; deba@2040: } deba@2040: bnode_matching[graph->bNode(it)] = it; deba@2040: } deba@2040: } deba@2040: return false; deba@2040: } deba@2040: deba@2040: /// \brief An augmenting phase of the Hopcroft-Karp algorithm deba@2040: /// deba@2040: /// It runs an augmenting phase of the Hopcroft-Karp deba@2040: /// algorithm. The phase finds maximum count of edge disjoint deba@2040: /// augmenting paths and augments on these paths. The algorithm deba@2040: /// consists at most of \f$ O(\sqrt{n}) \f$ phase and one phase is deba@2040: /// \f$ O(e) \f$ long. deba@2040: bool augment() { deba@2040: deba@2040: typename Graph::template ANodeMap areached(*graph, false); deba@2040: typename Graph::template BNodeMap breached(*graph, false); deba@2040: deba@2040: typename Graph::template BNodeMap bpred(*graph, INVALID); deba@2040: deba@2040: std::vector queue, bqueue; deba@2040: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2040: if (anode_matching[it] == INVALID) { deba@2040: queue.push_back(it); deba@2040: areached[it] = true; deba@2040: } deba@2040: } deba@2040: deba@2040: bool success = false; deba@2040: deba@2040: while (!success && !queue.empty()) { deba@2040: std::vector newqueue; deba@2040: for (int i = 0; i < (int)queue.size(); ++i) { deba@2040: Node anode = queue[i]; deba@2040: for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) { deba@2040: Node bnode = graph->bNode(jt); deba@2040: if (breached[bnode]) continue; deba@2040: breached[bnode] = true; deba@2040: bpred[bnode] = jt; deba@2040: if (bnode_matching[bnode] == INVALID) { deba@2040: bqueue.push_back(bnode); deba@2040: success = true; deba@2040: } else { deba@2040: Node newanode = graph->aNode(bnode_matching[bnode]); deba@2040: if (!areached[newanode]) { deba@2040: areached[newanode] = true; deba@2040: newqueue.push_back(newanode); deba@2040: } deba@2040: } deba@2040: } deba@2040: } deba@2040: queue.swap(newqueue); deba@2040: } deba@2040: deba@2040: if (success) { deba@2040: deba@2040: typename Graph::template ANodeMap aused(*graph, false); deba@2040: deba@2040: for (int i = 0; i < (int)bqueue.size(); ++i) { deba@2040: Node bnode = bqueue[i]; deba@2040: deba@2040: bool used = false; deba@2040: deba@2040: while (bnode != INVALID) { deba@2040: UEdge uedge = bpred[bnode]; deba@2040: Node anode = graph->aNode(uedge); deba@2040: deba@2040: if (aused[anode]) { deba@2040: used = true; deba@2040: break; deba@2040: } deba@2040: deba@2040: bnode = anode_matching[anode] != INVALID ? deba@2040: graph->bNode(anode_matching[anode]) : INVALID; deba@2040: deba@2040: } deba@2040: deba@2040: if (used) continue; deba@2040: deba@2040: bnode = bqueue[i]; deba@2040: while (bnode != INVALID) { deba@2040: UEdge uedge = bpred[bnode]; deba@2040: Node anode = graph->aNode(uedge); deba@2040: deba@2040: bnode_matching[bnode] = uedge; deba@2040: deba@2040: bnode = anode_matching[anode] != INVALID ? deba@2040: graph->bNode(anode_matching[anode]) : INVALID; deba@2040: deba@2040: anode_matching[anode] = uedge; deba@2040: deba@2040: aused[anode] = true; deba@2040: } deba@2051: ++matching_size; deba@2040: deba@2040: } deba@2040: } deba@2040: return success; deba@2040: } deba@2040: deba@2040: /// \brief An augmenting phase of the Ford-Fulkerson algorithm deba@2040: /// deba@2040: /// It runs an augmenting phase of the Ford-Fulkerson deba@2040: /// algorithm. The phase finds only one augmenting path and deba@2040: /// augments only on this paths. The algorithm consists at most deba@2040: /// of \f$ O(n) \f$ simple phase and one phase is at most deba@2040: /// \f$ O(e) \f$ long. deba@2040: bool simpleAugment() { deba@2040: deba@2040: typename Graph::template ANodeMap areached(*graph, false); deba@2040: typename Graph::template BNodeMap breached(*graph, false); deba@2040: deba@2040: typename Graph::template BNodeMap bpred(*graph, INVALID); deba@2040: deba@2040: std::vector queue; deba@2040: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2040: if (anode_matching[it] == INVALID) { deba@2040: queue.push_back(it); deba@2040: areached[it] = true; deba@2040: } deba@2040: } deba@2040: deba@2040: while (!queue.empty()) { deba@2040: std::vector newqueue; deba@2040: for (int i = 0; i < (int)queue.size(); ++i) { deba@2040: Node anode = queue[i]; deba@2040: for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) { deba@2040: Node bnode = graph->bNode(jt); deba@2040: if (breached[bnode]) continue; deba@2040: breached[bnode] = true; deba@2040: bpred[bnode] = jt; deba@2040: if (bnode_matching[bnode] == INVALID) { deba@2040: while (bnode != INVALID) { deba@2040: UEdge uedge = bpred[bnode]; deba@2040: anode = graph->aNode(uedge); deba@2040: deba@2040: bnode_matching[bnode] = uedge; deba@2040: deba@2040: bnode = anode_matching[anode] != INVALID ? deba@2040: graph->bNode(anode_matching[anode]) : INVALID; deba@2040: deba@2040: anode_matching[anode] = uedge; deba@2040: deba@2040: } deba@2051: ++matching_size; deba@2040: return true; deba@2040: } else { deba@2040: Node newanode = graph->aNode(bnode_matching[bnode]); deba@2040: if (!areached[newanode]) { deba@2040: areached[newanode] = true; deba@2040: newqueue.push_back(newanode); deba@2040: } deba@2040: } deba@2040: } deba@2040: } deba@2040: queue.swap(newqueue); deba@2040: } deba@2040: deba@2040: return false; deba@2040: } deba@2040: deba@2040: /// \brief Starts the algorithm. deba@2040: /// deba@2040: /// Starts the algorithm. It runs augmenting phases until the optimal deba@2040: /// solution reached. deba@2040: void start() { deba@2040: while (augment()) {} deba@2040: } deba@2040: deba@2040: /// \brief Runs the algorithm. deba@2040: /// deba@2040: /// It just initalize the algorithm and then start it. deba@2040: void run() { deba@2058: greedyInit(); deba@2040: start(); deba@2040: } deba@2040: deba@2040: /// @} deba@2040: deba@2040: /// \name Query Functions deba@2040: /// The result of the %Matching algorithm can be obtained using these deba@2040: /// functions.\n deba@2040: /// Before the use of these functions, deba@2040: /// either run() or start() must be called. deba@2040: deba@2040: ///@{ deba@2040: deba@2040: /// \brief Returns an minimum covering of the nodes. deba@2040: /// deba@2040: /// The minimum covering set problem is the dual solution of the deba@2040: /// maximum bipartite matching. It provides an solution for this deba@2040: /// problem what is proof of the optimality of the matching. deba@2040: /// \return The size of the cover set. deba@2040: template deba@2058: int coverSet(CoverMap& covering) const { deba@2040: deba@2040: typename Graph::template ANodeMap areached(*graph, false); deba@2040: typename Graph::template BNodeMap breached(*graph, false); deba@2040: deba@2040: std::vector queue; deba@2040: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2040: if (anode_matching[it] == INVALID) { deba@2040: queue.push_back(it); deba@2040: } deba@2040: } deba@2040: deba@2040: while (!queue.empty()) { deba@2040: std::vector newqueue; deba@2040: for (int i = 0; i < (int)queue.size(); ++i) { deba@2040: Node anode = queue[i]; deba@2040: for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) { deba@2040: Node bnode = graph->bNode(jt); deba@2040: if (breached[bnode]) continue; deba@2040: breached[bnode] = true; deba@2040: if (bnode_matching[bnode] != INVALID) { deba@2040: Node newanode = graph->aNode(bnode_matching[bnode]); deba@2040: if (!areached[newanode]) { deba@2040: areached[newanode] = true; deba@2040: newqueue.push_back(newanode); deba@2040: } deba@2040: } deba@2040: } deba@2040: } deba@2040: queue.swap(newqueue); deba@2040: } deba@2040: deba@2040: int size = 0; deba@2040: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2040: covering[it] = !areached[it] && anode_matching[it] != INVALID; deba@2040: if (!areached[it] && anode_matching[it] != INVALID) { deba@2040: ++size; deba@2040: } deba@2040: } deba@2040: for (BNodeIt it(*graph); it != INVALID; ++it) { deba@2040: covering[it] = breached[it]; deba@2040: if (breached[it]) { deba@2040: ++size; deba@2040: } deba@2040: } deba@2040: return size; deba@2040: } deba@2040: deba@2040: /// \brief Set true all matching uedge in the map. deba@2040: /// deba@2040: /// Set true all matching uedge in the map. It does not change the deba@2040: /// value mapped to the other uedges. deba@2040: /// \return The number of the matching edges. deba@2040: template deba@2058: int quickMatching(MatchingMap& matching) const { deba@2040: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2040: if (anode_matching[it] != INVALID) { deba@2040: matching[anode_matching[it]] = true; deba@2040: } deba@2040: } deba@2051: return matching_size; deba@2040: } deba@2040: deba@2040: /// \brief Set true all matching uedge in the map and the others to false. deba@2040: /// deba@2040: /// Set true all matching uedge in the map and the others to false. deba@2040: /// \return The number of the matching edges. deba@2040: template deba@2058: int matching(MatchingMap& matching) const { deba@2040: for (UEdgeIt it(*graph); it != INVALID; ++it) { deba@2040: matching[it] = it == anode_matching[graph->aNode(it)]; deba@2040: } deba@2051: return matching_size; deba@2040: } deba@2040: deba@2040: deba@2040: /// \brief Return true if the given uedge is in the matching. deba@2040: /// deba@2040: /// It returns true if the given uedge is in the matching. deba@2058: bool matchingEdge(const UEdge& edge) const { deba@2040: return anode_matching[graph->aNode(edge)] == edge; deba@2040: } deba@2040: deba@2040: /// \brief Returns the matching edge from the node. deba@2040: /// deba@2040: /// Returns the matching edge from the node. If there is not such deba@2040: /// edge it gives back \c INVALID. deba@2058: UEdge matchingEdge(const Node& node) const { deba@2040: if (graph->aNode(node)) { deba@2040: return anode_matching[node]; deba@2040: } else { deba@2040: return bnode_matching[node]; deba@2040: } deba@2040: } deba@2040: deba@2040: /// \brief Gives back the number of the matching edges. deba@2040: /// deba@2040: /// Gives back the number of the matching edges. deba@2051: int matchingSize() const { deba@2051: return matching_size; deba@2040: } deba@2040: deba@2040: /// @} deba@2040: deba@2040: private: deba@2040: deba@2040: ANodeMatchingMap anode_matching; deba@2040: BNodeMatchingMap bnode_matching; deba@2040: const Graph *graph; deba@2040: deba@2051: int matching_size; deba@2051: deba@2051: }; deba@2051: deba@2058: /// \ingroup matching deba@2058: /// deba@2058: /// \brief Maximum cardinality bipartite matching deba@2058: /// deba@2058: /// This function calculates the maximum cardinality matching deba@2058: /// in a bipartite graph. It gives back the matching in an undirected deba@2058: /// edge map. deba@2058: /// deba@2058: /// \param graph The bipartite graph. deba@2058: /// \retval matching The undirected edge map which will be set to deba@2058: /// the matching. deba@2058: /// \return The size of the matching. deba@2058: template deba@2058: int maxBipartiteMatching(const BpUGraph& graph, MatchingMap& matching) { deba@2058: MaxBipartiteMatching bpmatching(graph); deba@2058: bpmatching.run(); deba@2058: bpmatching.matching(matching); deba@2058: return bpmatching.matchingSize(); deba@2058: } deba@2058: deba@2051: /// \brief Default traits class for weighted bipartite matching algoritms. deba@2051: /// deba@2051: /// Default traits class for weighted bipartite matching algoritms. deba@2051: /// \param _BpUGraph The bipartite undirected graph type. deba@2051: /// \param _WeightMap Type of weight map. deba@2051: template deba@2051: struct WeightedBipartiteMatchingDefaultTraits { deba@2051: /// \brief The type of the weight of the undirected edges. deba@2051: typedef typename _WeightMap::Value Value; deba@2051: deba@2051: /// The undirected bipartite graph type the algorithm runs on. deba@2051: typedef _BpUGraph BpUGraph; deba@2051: deba@2051: /// The map of the edges weights deba@2051: typedef _WeightMap WeightMap; deba@2051: deba@2051: /// \brief The cross reference type used by heap. deba@2051: /// deba@2051: /// The cross reference type used by heap. deba@2051: /// Usually it is \c Graph::NodeMap. deba@2051: typedef typename BpUGraph::template NodeMap HeapCrossRef; deba@2051: deba@2051: /// \brief Instantiates a HeapCrossRef. deba@2051: /// deba@2051: /// This function instantiates a \ref HeapCrossRef. deba@2051: /// \param graph is the graph, to which we would like to define the deba@2051: /// HeapCrossRef. deba@2051: static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) { deba@2051: return new HeapCrossRef(graph); deba@2051: } deba@2051: deba@2051: /// \brief The heap type used by weighted matching algorithms. deba@2051: /// deba@2051: /// The heap type used by weighted matching algorithms. It should deba@2051: /// minimize the priorities and the heap's key type is the graph's deba@2051: /// anode graph's node. deba@2051: /// deba@2051: /// \sa BinHeap deba@2051: typedef BinHeap Heap; deba@2051: deba@2051: /// \brief Instantiates a Heap. deba@2051: /// deba@2051: /// This function instantiates a \ref Heap. deba@2051: /// \param crossref The cross reference of the heap. deba@2051: static Heap *createHeap(HeapCrossRef& crossref) { deba@2051: return new Heap(crossref); deba@2051: } deba@2051: deba@2051: }; deba@2051: deba@2051: deba@2051: /// \ingroup matching deba@2051: /// deba@2051: /// \brief Bipartite Max Weighted Matching algorithm deba@2051: /// deba@2051: /// This class implements the bipartite Max Weighted Matching deba@2051: /// algorithm. It uses the successive shortest path algorithm to deba@2051: /// calculate the maximum weighted matching in the bipartite deba@2051: /// graph. The algorithm can be used also to calculate the maximum deba@2051: /// cardinality maximum weighted matching. The time complexity deba@2051: /// of the algorithm is \f$ O(ne\log(n)) \f$ with the default binary deba@2051: /// heap implementation but this can be improved to deba@2051: /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps. deba@2051: /// deba@2051: /// The algorithm also provides a potential function on the nodes deba@2051: /// which a dual solution of the matching algorithm and it can be deba@2051: /// used to proof the optimality of the given pimal solution. deba@2051: #ifdef DOXYGEN deba@2051: template deba@2051: #else deba@2051: template , deba@2051: typename _Traits = WeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> > deba@2051: #endif deba@2051: class MaxWeightedBipartiteMatching { deba@2051: public: deba@2051: deba@2051: typedef _Traits Traits; deba@2051: typedef typename Traits::BpUGraph BpUGraph; deba@2051: typedef typename Traits::WeightMap WeightMap; deba@2051: typedef typename Traits::Value Value; deba@2051: deba@2051: protected: deba@2051: deba@2051: typedef typename Traits::HeapCrossRef HeapCrossRef; deba@2051: typedef typename Traits::Heap Heap; deba@2051: deba@2051: deba@2051: typedef typename BpUGraph::Node Node; deba@2051: typedef typename BpUGraph::ANodeIt ANodeIt; deba@2051: typedef typename BpUGraph::BNodeIt BNodeIt; deba@2051: typedef typename BpUGraph::UEdge UEdge; deba@2051: typedef typename BpUGraph::UEdgeIt UEdgeIt; deba@2051: typedef typename BpUGraph::IncEdgeIt IncEdgeIt; deba@2051: deba@2051: typedef typename BpUGraph::template ANodeMap ANodeMatchingMap; deba@2051: typedef typename BpUGraph::template BNodeMap BNodeMatchingMap; deba@2051: deba@2051: typedef typename BpUGraph::template ANodeMap ANodePotentialMap; deba@2051: typedef typename BpUGraph::template BNodeMap BNodePotentialMap; deba@2051: deba@2051: deba@2051: public: deba@2051: deba@2051: /// \brief \ref Exception for uninitialized parameters. deba@2051: /// deba@2051: /// This error represents problems in the initialization deba@2051: /// of the parameters of the algorithms. deba@2051: class UninitializedParameter : public lemon::UninitializedParameter { deba@2051: public: alpar@2151: virtual const char* what() const throw() { deba@2051: return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter"; deba@2051: } deba@2051: }; deba@2051: deba@2051: ///\name Named template parameters deba@2051: deba@2051: ///@{ deba@2051: deba@2051: template deba@2051: struct DefHeapTraits : public Traits { deba@2051: typedef CR HeapCrossRef; deba@2051: typedef H Heap; deba@2051: static HeapCrossRef *createHeapCrossRef(const BpUGraph &) { deba@2051: throw UninitializedParameter(); deba@2051: } deba@2051: static Heap *createHeap(HeapCrossRef &) { deba@2051: throw UninitializedParameter(); deba@2051: } deba@2051: }; deba@2051: deba@2051: /// \brief \ref named-templ-param "Named parameter" for setting heap deba@2051: /// and cross reference type deba@2051: /// deba@2051: /// \ref named-templ-param "Named parameter" for setting heap and cross deba@2051: /// reference type deba@2051: template > deba@2051: struct DefHeap deba@2051: : public MaxWeightedBipartiteMatching > { deba@2051: typedef MaxWeightedBipartiteMatching > Create; deba@2051: }; deba@2051: deba@2051: template deba@2051: struct DefStandardHeapTraits : public Traits { deba@2051: typedef CR HeapCrossRef; deba@2051: typedef H Heap; deba@2051: static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) { deba@2051: return new HeapCrossRef(graph); deba@2051: } deba@2051: static Heap *createHeap(HeapCrossRef &crossref) { deba@2051: return new Heap(crossref); deba@2051: } deba@2051: }; deba@2051: deba@2051: /// \brief \ref named-templ-param "Named parameter" for setting heap and deba@2051: /// cross reference type with automatic allocation deba@2051: /// deba@2051: /// \ref named-templ-param "Named parameter" for setting heap and cross deba@2051: /// reference type. It can allocate the heap and the cross reference deba@2051: /// object if the cross reference's constructor waits for the graph as deba@2051: /// parameter and the heap's constructor waits for the cross reference. deba@2051: template > deba@2051: struct DefStandardHeap deba@2051: : public MaxWeightedBipartiteMatching > { deba@2051: typedef MaxWeightedBipartiteMatching > deba@2051: Create; deba@2051: }; deba@2051: deba@2051: ///@} deba@2051: deba@2051: deba@2051: /// \brief Constructor. deba@2051: /// deba@2051: /// Constructor of the algorithm. deba@2051: MaxWeightedBipartiteMatching(const BpUGraph& _graph, deba@2051: const WeightMap& _weight) deba@2051: : graph(&_graph), weight(&_weight), deba@2051: anode_matching(_graph), bnode_matching(_graph), deba@2051: anode_potential(_graph), bnode_potential(_graph), deba@2051: _heap_cross_ref(0), local_heap_cross_ref(false), deba@2051: _heap(0), local_heap(0) {} deba@2051: deba@2051: /// \brief Destructor. deba@2051: /// deba@2051: /// Destructor of the algorithm. deba@2051: ~MaxWeightedBipartiteMatching() { deba@2051: destroyStructures(); deba@2051: } deba@2051: deba@2051: /// \brief Sets the heap and the cross reference used by algorithm. deba@2051: /// deba@2051: /// Sets the heap and the cross reference used by algorithm. deba@2051: /// If you don't use this function before calling \ref run(), deba@2051: /// it will allocate one. The destuctor deallocates this deba@2051: /// automatically allocated map, of course. deba@2051: /// \return \c (*this) deba@2051: MaxWeightedBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) { deba@2051: if(local_heap_cross_ref) { deba@2051: delete _heap_cross_ref; deba@2051: local_heap_cross_ref = false; deba@2051: } deba@2051: _heap_cross_ref = &crossRef; deba@2051: if(local_heap) { deba@2051: delete _heap; deba@2051: local_heap = false; deba@2051: } deba@2051: _heap = &heap; deba@2051: return *this; deba@2051: } deba@2051: deba@2051: /// \name Execution control deba@2051: /// The simplest way to execute the algorithm is to use deba@2051: /// one of the member functions called \c run(). deba@2051: /// \n deba@2051: /// If you need more control on the execution, deba@2051: /// first you must call \ref init() or one alternative for it. deba@2051: /// Finally \ref start() will perform the matching computation or deba@2051: /// with step-by-step execution you can augment the solution. deba@2051: deba@2051: /// @{ deba@2051: deba@2051: /// \brief Initalize the data structures. deba@2051: /// deba@2051: /// It initalizes the data structures and creates an empty matching. deba@2051: void init() { deba@2051: initStructures(); deba@2051: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2051: anode_matching[it] = INVALID; deba@2051: anode_potential[it] = 0; deba@2051: } deba@2051: for (BNodeIt it(*graph); it != INVALID; ++it) { deba@2051: bnode_matching[it] = INVALID; deba@2051: bnode_potential[it] = 0; deba@2051: for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) { deba@2058: if ((*weight)[jt] > bnode_potential[it]) { deba@2058: bnode_potential[it] = (*weight)[jt]; deba@2051: } deba@2051: } deba@2051: } deba@2051: matching_value = 0; deba@2051: matching_size = 0; deba@2051: } deba@2051: deba@2051: deba@2051: /// \brief An augmenting phase of the weighted matching algorithm deba@2051: /// deba@2051: /// It runs an augmenting phase of the weighted matching deba@2051: /// algorithm. The phase finds the best augmenting path and deba@2051: /// augments only on this paths. deba@2051: /// deba@2051: /// The algorithm consists at most deba@2051: /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ deba@2051: /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long deba@2051: /// with binary heap. deba@2051: /// \param decrease If the given parameter true the matching value deba@2051: /// can be decreased in the augmenting phase. If we would like deba@2051: /// to calculate the maximum cardinality maximum weighted matching deba@2051: /// then we should let the algorithm to decrease the matching deba@2051: /// value in order to increase the number of the matching edges. deba@2051: bool augment(bool decrease = false) { deba@2051: deba@2051: typename BpUGraph::template BNodeMap bdist(*graph); deba@2051: typename BpUGraph::template BNodeMap bpred(*graph, INVALID); deba@2051: deba@2051: Node bestNode = INVALID; deba@2051: Value bestValue = 0; deba@2051: deba@2051: _heap->clear(); deba@2051: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2051: (*_heap_cross_ref)[it] = Heap::PRE_HEAP; deba@2051: } deba@2051: deba@2051: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2051: if (anode_matching[it] == INVALID) { deba@2051: _heap->push(it, 0); deba@2051: } deba@2051: } deba@2051: deba@2051: Value bdistMax = 0; deba@2051: while (!_heap->empty()) { deba@2051: Node anode = _heap->top(); deba@2051: Value avalue = _heap->prio(); deba@2051: _heap->pop(); deba@2051: for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) { deba@2051: if (jt == anode_matching[anode]) continue; deba@2051: Node bnode = graph->bNode(jt); deba@2058: Value bvalue = avalue - (*weight)[jt] + deba@2058: anode_potential[anode] + bnode_potential[bnode]; deba@2051: if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) { deba@2051: bdist[bnode] = bvalue; deba@2051: bpred[bnode] = jt; deba@2051: } deba@2051: if (bvalue > bdistMax) { deba@2051: bdistMax = bvalue; deba@2051: } deba@2051: if (bnode_matching[bnode] != INVALID) { deba@2051: Node newanode = graph->aNode(bnode_matching[bnode]); deba@2051: switch (_heap->state(newanode)) { deba@2051: case Heap::PRE_HEAP: deba@2051: _heap->push(newanode, bvalue); deba@2051: break; deba@2051: case Heap::IN_HEAP: deba@2051: if (bvalue < (*_heap)[newanode]) { deba@2051: _heap->decrease(newanode, bvalue); deba@2051: } deba@2051: break; deba@2051: case Heap::POST_HEAP: deba@2051: break; deba@2051: } deba@2051: } else { deba@2051: if (bestNode == INVALID || deba@2058: bnode_potential[bnode] - bvalue > bestValue) { deba@2058: bestValue = bnode_potential[bnode] - bvalue; deba@2051: bestNode = bnode; deba@2051: } deba@2051: } deba@2051: } deba@2051: } deba@2051: deba@2051: if (bestNode == INVALID || (!decrease && bestValue < 0)) { deba@2051: return false; deba@2051: } deba@2051: deba@2051: matching_value += bestValue; deba@2051: ++matching_size; deba@2051: deba@2051: for (BNodeIt it(*graph); it != INVALID; ++it) { deba@2051: if (bpred[it] != INVALID) { deba@2058: bnode_potential[it] -= bdist[it]; deba@2051: } else { deba@2058: bnode_potential[it] -= bdistMax; deba@2051: } deba@2051: } deba@2051: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2051: if (anode_matching[it] != INVALID) { deba@2051: Node bnode = graph->bNode(anode_matching[it]); deba@2051: if (bpred[bnode] != INVALID) { deba@2051: anode_potential[it] += bdist[bnode]; deba@2051: } else { deba@2051: anode_potential[it] += bdistMax; deba@2051: } deba@2051: } deba@2051: } deba@2051: deba@2051: while (bestNode != INVALID) { deba@2051: UEdge uedge = bpred[bestNode]; deba@2051: Node anode = graph->aNode(uedge); deba@2051: deba@2051: bnode_matching[bestNode] = uedge; deba@2051: if (anode_matching[anode] != INVALID) { deba@2051: bestNode = graph->bNode(anode_matching[anode]); deba@2051: } else { deba@2051: bestNode = INVALID; deba@2051: } deba@2051: anode_matching[anode] = uedge; deba@2051: } deba@2051: deba@2051: deba@2051: return true; deba@2051: } deba@2051: deba@2051: /// \brief Starts the algorithm. deba@2051: /// deba@2051: /// Starts the algorithm. It runs augmenting phases until the deba@2051: /// optimal solution reached. deba@2051: /// deba@2051: /// \param maxCardinality If the given value is true it will deba@2051: /// calculate the maximum cardinality maximum matching instead of deba@2051: /// the maximum matching. deba@2051: void start(bool maxCardinality = false) { deba@2051: while (augment(maxCardinality)) {} deba@2051: } deba@2051: deba@2051: /// \brief Runs the algorithm. deba@2051: /// deba@2051: /// It just initalize the algorithm and then start it. deba@2051: /// deba@2051: /// \param maxCardinality If the given value is true it will deba@2051: /// calculate the maximum cardinality maximum matching instead of deba@2051: /// the maximum matching. deba@2051: void run(bool maxCardinality = false) { deba@2051: init(); deba@2051: start(maxCardinality); deba@2051: } deba@2051: deba@2051: /// @} deba@2051: deba@2051: /// \name Query Functions deba@2051: /// The result of the %Matching algorithm can be obtained using these deba@2051: /// functions.\n deba@2051: /// Before the use of these functions, deba@2051: /// either run() or start() must be called. deba@2051: deba@2051: ///@{ deba@2051: deba@2051: /// \brief Gives back the potential in the NodeMap deba@2051: /// deba@2058: /// Gives back the potential in the NodeMap. The matching is optimal deba@2058: /// with the current number of edges if \f$ \pi(a) + \pi(b) - w(ab) = 0 \f$ deba@2058: /// for each matching edges and \f$ \pi(a) + \pi(b) - w(ab) \ge 0 \f$ deba@2058: /// for each edges. deba@2051: template deba@2058: void potential(PotentialMap& potential) const { deba@2051: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2051: potential[it] = anode_potential[it]; deba@2051: } deba@2051: for (BNodeIt it(*graph); it != INVALID; ++it) { deba@2051: potential[it] = bnode_potential[it]; deba@2051: } deba@2051: } deba@2051: deba@2051: /// \brief Set true all matching uedge in the map. deba@2051: /// deba@2051: /// Set true all matching uedge in the map. It does not change the deba@2051: /// value mapped to the other uedges. deba@2051: /// \return The number of the matching edges. deba@2051: template deba@2058: int quickMatching(MatchingMap& matching) const { deba@2051: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2051: if (anode_matching[it] != INVALID) { deba@2051: matching[anode_matching[it]] = true; deba@2051: } deba@2051: } deba@2051: return matching_size; deba@2051: } deba@2051: deba@2051: /// \brief Set true all matching uedge in the map and the others to false. deba@2051: /// deba@2051: /// Set true all matching uedge in the map and the others to false. deba@2051: /// \return The number of the matching edges. deba@2051: template deba@2058: int matching(MatchingMap& matching) const { deba@2051: for (UEdgeIt it(*graph); it != INVALID; ++it) { deba@2051: matching[it] = it == anode_matching[graph->aNode(it)]; deba@2051: } deba@2051: return matching_size; deba@2051: } deba@2051: deba@2051: deba@2051: /// \brief Return true if the given uedge is in the matching. deba@2051: /// deba@2051: /// It returns true if the given uedge is in the matching. deba@2058: bool matchingEdge(const UEdge& edge) const { deba@2051: return anode_matching[graph->aNode(edge)] == edge; deba@2051: } deba@2051: deba@2051: /// \brief Returns the matching edge from the node. deba@2051: /// deba@2051: /// Returns the matching edge from the node. If there is not such deba@2051: /// edge it gives back \c INVALID. deba@2058: UEdge matchingEdge(const Node& node) const { deba@2051: if (graph->aNode(node)) { deba@2051: return anode_matching[node]; deba@2051: } else { deba@2051: return bnode_matching[node]; deba@2051: } deba@2051: } deba@2051: deba@2051: /// \brief Gives back the sum of weights of the matching edges. deba@2051: /// deba@2051: /// Gives back the sum of weights of the matching edges. deba@2051: Value matchingValue() const { deba@2051: return matching_value; deba@2051: } deba@2051: deba@2051: /// \brief Gives back the number of the matching edges. deba@2051: /// deba@2051: /// Gives back the number of the matching edges. deba@2051: int matchingSize() const { deba@2051: return matching_size; deba@2051: } deba@2051: deba@2051: /// @} deba@2051: deba@2051: private: deba@2051: deba@2051: void initStructures() { deba@2051: if (!_heap_cross_ref) { deba@2051: local_heap_cross_ref = true; deba@2051: _heap_cross_ref = Traits::createHeapCrossRef(*graph); deba@2051: } deba@2051: if (!_heap) { deba@2051: local_heap = true; deba@2051: _heap = Traits::createHeap(*_heap_cross_ref); deba@2051: } deba@2051: } deba@2051: deba@2051: void destroyStructures() { deba@2051: if (local_heap_cross_ref) delete _heap_cross_ref; deba@2051: if (local_heap) delete _heap; deba@2051: } deba@2051: deba@2051: deba@2051: private: deba@2051: deba@2051: const BpUGraph *graph; deba@2051: const WeightMap* weight; deba@2051: deba@2051: ANodeMatchingMap anode_matching; deba@2051: BNodeMatchingMap bnode_matching; deba@2051: deba@2051: ANodePotentialMap anode_potential; deba@2051: BNodePotentialMap bnode_potential; deba@2051: deba@2051: Value matching_value; deba@2051: int matching_size; deba@2051: deba@2051: HeapCrossRef *_heap_cross_ref; deba@2051: bool local_heap_cross_ref; deba@2051: deba@2051: Heap *_heap; deba@2051: bool local_heap; deba@2051: deba@2051: }; deba@2051: deba@2058: /// \ingroup matching deba@2058: /// deba@2058: /// \brief Maximum weighted bipartite matching deba@2058: /// deba@2058: /// This function calculates the maximum weighted matching deba@2058: /// in a bipartite graph. It gives back the matching in an undirected deba@2058: /// edge map. deba@2058: /// deba@2058: /// \param graph The bipartite graph. deba@2058: /// \param weight The undirected edge map which contains the weights. deba@2058: /// \retval matching The undirected edge map which will be set to deba@2058: /// the matching. deba@2058: /// \return The value of the matching. deba@2058: template deba@2058: typename WeightMap::Value deba@2058: maxWeightedBipartiteMatching(const BpUGraph& graph, const WeightMap& weight, deba@2058: MatchingMap& matching) { deba@2058: MaxWeightedBipartiteMatching deba@2058: bpmatching(graph, weight); deba@2058: bpmatching.run(); deba@2058: bpmatching.matching(matching); deba@2058: return bpmatching.matchingValue(); deba@2058: } deba@2058: deba@2058: /// \ingroup matching deba@2058: /// deba@2058: /// \brief Maximum weighted maximum cardinality bipartite matching deba@2058: /// deba@2058: /// This function calculates the maximum weighted of the maximum cardinality deba@2058: /// matchings of a bipartite graph. It gives back the matching in an deba@2058: /// undirected edge map. deba@2058: /// deba@2058: /// \param graph The bipartite graph. deba@2058: /// \param weight The undirected edge map which contains the weights. deba@2058: /// \retval matching The undirected edge map which will be set to deba@2058: /// the matching. deba@2058: /// \return The value of the matching. deba@2058: template deba@2058: typename WeightMap::Value deba@2058: maxWeightedMaxBipartiteMatching(const BpUGraph& graph, deba@2058: const WeightMap& weight, deba@2058: MatchingMap& matching) { deba@2058: MaxWeightedBipartiteMatching deba@2058: bpmatching(graph, weight); deba@2058: bpmatching.run(true); deba@2058: bpmatching.matching(matching); deba@2058: return bpmatching.matchingValue(); deba@2058: } deba@2058: deba@2051: /// \brief Default traits class for minimum cost bipartite matching deba@2051: /// algoritms. deba@2051: /// deba@2051: /// Default traits class for minimum cost bipartite matching deba@2051: /// algoritms. deba@2051: /// deba@2051: /// \param _BpUGraph The bipartite undirected graph deba@2051: /// type. deba@2051: /// deba@2051: /// \param _CostMap Type of cost map. deba@2051: template deba@2051: struct MinCostMaxBipartiteMatchingDefaultTraits { deba@2051: /// \brief The type of the cost of the undirected edges. deba@2051: typedef typename _CostMap::Value Value; deba@2051: deba@2051: /// The undirected bipartite graph type the algorithm runs on. deba@2051: typedef _BpUGraph BpUGraph; deba@2051: deba@2051: /// The map of the edges costs deba@2051: typedef _CostMap CostMap; deba@2051: deba@2051: /// \brief The cross reference type used by heap. deba@2051: /// deba@2051: /// The cross reference type used by heap. deba@2051: /// Usually it is \c Graph::NodeMap. deba@2051: typedef typename BpUGraph::template NodeMap HeapCrossRef; deba@2051: deba@2051: /// \brief Instantiates a HeapCrossRef. deba@2051: /// deba@2051: /// This function instantiates a \ref HeapCrossRef. deba@2051: /// \param graph is the graph, to which we would like to define the deba@2051: /// HeapCrossRef. deba@2051: static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) { deba@2051: return new HeapCrossRef(graph); deba@2051: } deba@2051: deba@2051: /// \brief The heap type used by costed matching algorithms. deba@2051: /// deba@2051: /// The heap type used by costed matching algorithms. It should deba@2051: /// minimize the priorities and the heap's key type is the graph's deba@2051: /// anode graph's node. deba@2051: /// deba@2051: /// \sa BinHeap deba@2051: typedef BinHeap Heap; deba@2051: deba@2051: /// \brief Instantiates a Heap. deba@2051: /// deba@2051: /// This function instantiates a \ref Heap. deba@2051: /// \param crossref The cross reference of the heap. deba@2051: static Heap *createHeap(HeapCrossRef& crossref) { deba@2051: return new Heap(crossref); deba@2051: } deba@2051: deba@2051: }; deba@2051: deba@2051: deba@2051: /// \ingroup matching deba@2051: /// deba@2051: /// \brief Bipartite Min Cost Matching algorithm deba@2051: /// deba@2051: /// This class implements the bipartite Min Cost Matching algorithm. deba@2051: /// It uses the successive shortest path algorithm to calculate the deba@2051: /// minimum cost maximum matching in the bipartite graph. The time deba@2051: /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the deba@2051: /// default binary heap implementation but this can be improved to deba@2051: /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps. deba@2051: /// deba@2051: /// The algorithm also provides a potential function on the nodes deba@2051: /// which a dual solution of the matching algorithm and it can be deba@2051: /// used to proof the optimality of the given pimal solution. deba@2051: #ifdef DOXYGEN deba@2051: template deba@2051: #else deba@2051: template , deba@2051: typename _Traits = MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> > deba@2051: #endif deba@2051: class MinCostMaxBipartiteMatching { deba@2051: public: deba@2051: deba@2051: typedef _Traits Traits; deba@2051: typedef typename Traits::BpUGraph BpUGraph; deba@2051: typedef typename Traits::CostMap CostMap; deba@2051: typedef typename Traits::Value Value; deba@2051: deba@2051: protected: deba@2051: deba@2051: typedef typename Traits::HeapCrossRef HeapCrossRef; deba@2051: typedef typename Traits::Heap Heap; deba@2051: deba@2051: deba@2051: typedef typename BpUGraph::Node Node; deba@2051: typedef typename BpUGraph::ANodeIt ANodeIt; deba@2051: typedef typename BpUGraph::BNodeIt BNodeIt; deba@2051: typedef typename BpUGraph::UEdge UEdge; deba@2051: typedef typename BpUGraph::UEdgeIt UEdgeIt; deba@2051: typedef typename BpUGraph::IncEdgeIt IncEdgeIt; deba@2051: deba@2051: typedef typename BpUGraph::template ANodeMap ANodeMatchingMap; deba@2051: typedef typename BpUGraph::template BNodeMap BNodeMatchingMap; deba@2051: deba@2051: typedef typename BpUGraph::template ANodeMap ANodePotentialMap; deba@2051: typedef typename BpUGraph::template BNodeMap BNodePotentialMap; deba@2051: deba@2051: deba@2051: public: deba@2051: deba@2051: /// \brief \ref Exception for uninitialized parameters. deba@2051: /// deba@2051: /// This error represents problems in the initialization deba@2051: /// of the parameters of the algorithms. deba@2051: class UninitializedParameter : public lemon::UninitializedParameter { deba@2051: public: alpar@2151: virtual const char* what() const throw() { deba@2051: return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter"; deba@2051: } deba@2051: }; deba@2051: deba@2051: ///\name Named template parameters deba@2051: deba@2051: ///@{ deba@2051: deba@2051: template deba@2051: struct DefHeapTraits : public Traits { deba@2051: typedef CR HeapCrossRef; deba@2051: typedef H Heap; deba@2051: static HeapCrossRef *createHeapCrossRef(const BpUGraph &) { deba@2051: throw UninitializedParameter(); deba@2051: } deba@2051: static Heap *createHeap(HeapCrossRef &) { deba@2051: throw UninitializedParameter(); deba@2051: } deba@2051: }; deba@2051: deba@2051: /// \brief \ref named-templ-param "Named parameter" for setting heap deba@2051: /// and cross reference type deba@2051: /// deba@2051: /// \ref named-templ-param "Named parameter" for setting heap and cross deba@2051: /// reference type deba@2051: template > deba@2051: struct DefHeap deba@2051: : public MinCostMaxBipartiteMatching > { deba@2051: typedef MinCostMaxBipartiteMatching > Create; deba@2051: }; deba@2051: deba@2051: template deba@2051: struct DefStandardHeapTraits : public Traits { deba@2051: typedef CR HeapCrossRef; deba@2051: typedef H Heap; deba@2051: static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) { deba@2051: return new HeapCrossRef(graph); deba@2051: } deba@2051: static Heap *createHeap(HeapCrossRef &crossref) { deba@2051: return new Heap(crossref); deba@2051: } deba@2051: }; deba@2051: deba@2051: /// \brief \ref named-templ-param "Named parameter" for setting heap and deba@2051: /// cross reference type with automatic allocation deba@2051: /// deba@2051: /// \ref named-templ-param "Named parameter" for setting heap and cross deba@2051: /// reference type. It can allocate the heap and the cross reference deba@2051: /// object if the cross reference's constructor waits for the graph as deba@2051: /// parameter and the heap's constructor waits for the cross reference. deba@2051: template > deba@2051: struct DefStandardHeap deba@2051: : public MinCostMaxBipartiteMatching > { deba@2051: typedef MinCostMaxBipartiteMatching > deba@2051: Create; deba@2051: }; deba@2051: deba@2051: ///@} deba@2051: deba@2051: deba@2051: /// \brief Constructor. deba@2051: /// deba@2051: /// Constructor of the algorithm. deba@2051: MinCostMaxBipartiteMatching(const BpUGraph& _graph, deba@2051: const CostMap& _cost) deba@2051: : graph(&_graph), cost(&_cost), deba@2051: anode_matching(_graph), bnode_matching(_graph), deba@2051: anode_potential(_graph), bnode_potential(_graph), deba@2051: _heap_cross_ref(0), local_heap_cross_ref(false), deba@2051: _heap(0), local_heap(0) {} deba@2051: deba@2051: /// \brief Destructor. deba@2051: /// deba@2051: /// Destructor of the algorithm. deba@2051: ~MinCostMaxBipartiteMatching() { deba@2051: destroyStructures(); deba@2051: } deba@2051: deba@2051: /// \brief Sets the heap and the cross reference used by algorithm. deba@2051: /// deba@2051: /// Sets the heap and the cross reference used by algorithm. deba@2051: /// If you don't use this function before calling \ref run(), deba@2051: /// it will allocate one. The destuctor deallocates this deba@2051: /// automatically allocated map, of course. deba@2051: /// \return \c (*this) deba@2051: MinCostMaxBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) { deba@2051: if(local_heap_cross_ref) { deba@2051: delete _heap_cross_ref; deba@2051: local_heap_cross_ref = false; deba@2051: } deba@2051: _heap_cross_ref = &crossRef; deba@2051: if(local_heap) { deba@2051: delete _heap; deba@2051: local_heap = false; deba@2051: } deba@2051: _heap = &heap; deba@2051: return *this; deba@2051: } deba@2051: deba@2051: /// \name Execution control deba@2051: /// The simplest way to execute the algorithm is to use deba@2051: /// one of the member functions called \c run(). deba@2051: /// \n deba@2051: /// If you need more control on the execution, deba@2051: /// first you must call \ref init() or one alternative for it. deba@2051: /// Finally \ref start() will perform the matching computation or deba@2051: /// with step-by-step execution you can augment the solution. deba@2051: deba@2051: /// @{ deba@2051: deba@2051: /// \brief Initalize the data structures. deba@2051: /// deba@2051: /// It initalizes the data structures and creates an empty matching. deba@2051: void init() { deba@2051: initStructures(); deba@2051: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2051: anode_matching[it] = INVALID; deba@2051: anode_potential[it] = 0; deba@2051: } deba@2051: for (BNodeIt it(*graph); it != INVALID; ++it) { deba@2051: bnode_matching[it] = INVALID; deba@2051: bnode_potential[it] = 0; deba@2051: } deba@2051: matching_cost = 0; deba@2051: matching_size = 0; deba@2051: } deba@2051: deba@2051: deba@2051: /// \brief An augmenting phase of the costed matching algorithm deba@2051: /// deba@2051: /// It runs an augmenting phase of the matching algorithm. The deba@2051: /// phase finds the best augmenting path and augments only on this deba@2051: /// paths. deba@2051: /// deba@2051: /// The algorithm consists at most deba@2051: /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ deba@2051: /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long deba@2051: /// with binary heap. deba@2051: bool augment() { deba@2051: deba@2051: typename BpUGraph::template BNodeMap bdist(*graph); deba@2051: typename BpUGraph::template BNodeMap bpred(*graph, INVALID); deba@2051: deba@2051: Node bestNode = INVALID; deba@2051: Value bestValue = 0; deba@2051: deba@2051: _heap->clear(); deba@2051: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2051: (*_heap_cross_ref)[it] = Heap::PRE_HEAP; deba@2051: } deba@2051: deba@2051: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2051: if (anode_matching[it] == INVALID) { deba@2051: _heap->push(it, 0); deba@2051: } deba@2051: } deba@2136: Value bdistMax = 0; deba@2051: deba@2051: while (!_heap->empty()) { deba@2051: Node anode = _heap->top(); deba@2051: Value avalue = _heap->prio(); deba@2051: _heap->pop(); deba@2051: for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) { deba@2051: if (jt == anode_matching[anode]) continue; deba@2051: Node bnode = graph->bNode(jt); deba@2051: Value bvalue = avalue + (*cost)[jt] + deba@2051: anode_potential[anode] - bnode_potential[bnode]; deba@2051: if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) { deba@2051: bdist[bnode] = bvalue; deba@2051: bpred[bnode] = jt; deba@2051: } deba@2136: if (bvalue > bdistMax) { deba@2136: bdistMax = bvalue; deba@2136: } deba@2051: if (bnode_matching[bnode] != INVALID) { deba@2051: Node newanode = graph->aNode(bnode_matching[bnode]); deba@2051: switch (_heap->state(newanode)) { deba@2051: case Heap::PRE_HEAP: deba@2051: _heap->push(newanode, bvalue); deba@2051: break; deba@2051: case Heap::IN_HEAP: deba@2051: if (bvalue < (*_heap)[newanode]) { deba@2051: _heap->decrease(newanode, bvalue); deba@2051: } deba@2051: break; deba@2051: case Heap::POST_HEAP: deba@2051: break; deba@2051: } deba@2051: } else { deba@2051: if (bestNode == INVALID || deba@2051: bvalue + bnode_potential[bnode] < bestValue) { deba@2051: bestValue = bvalue + bnode_potential[bnode]; deba@2051: bestNode = bnode; deba@2051: } deba@2051: } deba@2051: } deba@2051: } deba@2051: deba@2051: if (bestNode == INVALID) { deba@2051: return false; deba@2051: } deba@2051: deba@2051: matching_cost += bestValue; deba@2051: ++matching_size; deba@2051: deba@2051: for (BNodeIt it(*graph); it != INVALID; ++it) { deba@2051: if (bpred[it] != INVALID) { deba@2051: bnode_potential[it] += bdist[it]; deba@2136: } else { deba@2136: bnode_potential[it] += bdistMax; deba@2051: } deba@2051: } deba@2051: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2051: if (anode_matching[it] != INVALID) { deba@2051: Node bnode = graph->bNode(anode_matching[it]); deba@2051: if (bpred[bnode] != INVALID) { deba@2051: anode_potential[it] += bdist[bnode]; deba@2136: } else { deba@2136: anode_potential[it] += bdistMax; deba@2051: } deba@2051: } deba@2051: } deba@2051: deba@2051: while (bestNode != INVALID) { deba@2051: UEdge uedge = bpred[bestNode]; deba@2051: Node anode = graph->aNode(uedge); deba@2051: deba@2051: bnode_matching[bestNode] = uedge; deba@2051: if (anode_matching[anode] != INVALID) { deba@2051: bestNode = graph->bNode(anode_matching[anode]); deba@2051: } else { deba@2051: bestNode = INVALID; deba@2051: } deba@2051: anode_matching[anode] = uedge; deba@2051: } deba@2051: deba@2051: deba@2051: return true; deba@2051: } deba@2051: deba@2051: /// \brief Starts the algorithm. deba@2051: /// deba@2051: /// Starts the algorithm. It runs augmenting phases until the deba@2051: /// optimal solution reached. deba@2051: void start() { deba@2051: while (augment()) {} deba@2051: } deba@2051: deba@2051: /// \brief Runs the algorithm. deba@2051: /// deba@2051: /// It just initalize the algorithm and then start it. deba@2051: void run() { deba@2051: init(); deba@2051: start(); deba@2051: } deba@2051: deba@2051: /// @} deba@2051: deba@2051: /// \name Query Functions deba@2051: /// The result of the %Matching algorithm can be obtained using these deba@2051: /// functions.\n deba@2051: /// Before the use of these functions, deba@2051: /// either run() or start() must be called. deba@2051: deba@2051: ///@{ deba@2051: deba@2051: /// \brief Gives back the potential in the NodeMap deba@2051: /// deba@2058: /// Gives back the potential in the NodeMap. The potential is optimal with deba@2058: /// the current number of edges if \f$ \pi(a) - \pi(b) + w(ab) = 0 \f$ for deba@2051: /// each matching edges and \f$ \pi(a) - \pi(b) + w(ab) \ge 0 \f$ deba@2051: /// for each edges. deba@2051: template deba@2058: void potential(PotentialMap& potential) const { deba@2051: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2051: potential[it] = anode_potential[it]; deba@2051: } deba@2051: for (BNodeIt it(*graph); it != INVALID; ++it) { deba@2051: potential[it] = bnode_potential[it]; deba@2051: } deba@2051: } deba@2051: deba@2051: /// \brief Set true all matching uedge in the map. deba@2051: /// deba@2051: /// Set true all matching uedge in the map. It does not change the deba@2051: /// value mapped to the other uedges. deba@2051: /// \return The number of the matching edges. deba@2051: template deba@2058: int quickMatching(MatchingMap& matching) const { deba@2051: for (ANodeIt it(*graph); it != INVALID; ++it) { deba@2051: if (anode_matching[it] != INVALID) { deba@2051: matching[anode_matching[it]] = true; deba@2051: } deba@2051: } deba@2051: return matching_size; deba@2051: } deba@2051: deba@2051: /// \brief Set true all matching uedge in the map and the others to false. deba@2051: /// deba@2051: /// Set true all matching uedge in the map and the others to false. deba@2051: /// \return The number of the matching edges. deba@2051: template deba@2058: int matching(MatchingMap& matching) const { deba@2051: for (UEdgeIt it(*graph); it != INVALID; ++it) { deba@2051: matching[it] = it == anode_matching[graph->aNode(it)]; deba@2051: } deba@2051: return matching_size; deba@2051: } deba@2051: deba@2051: deba@2051: /// \brief Return true if the given uedge is in the matching. deba@2051: /// deba@2051: /// It returns true if the given uedge is in the matching. deba@2058: bool matchingEdge(const UEdge& edge) const { deba@2051: return anode_matching[graph->aNode(edge)] == edge; deba@2051: } deba@2051: deba@2051: /// \brief Returns the matching edge from the node. deba@2051: /// deba@2051: /// Returns the matching edge from the node. If there is not such deba@2051: /// edge it gives back \c INVALID. deba@2058: UEdge matchingEdge(const Node& node) const { deba@2051: if (graph->aNode(node)) { deba@2051: return anode_matching[node]; deba@2051: } else { deba@2051: return bnode_matching[node]; deba@2051: } deba@2051: } deba@2051: deba@2051: /// \brief Gives back the sum of costs of the matching edges. deba@2051: /// deba@2051: /// Gives back the sum of costs of the matching edges. deba@2051: Value matchingCost() const { deba@2051: return matching_cost; deba@2051: } deba@2051: deba@2051: /// \brief Gives back the number of the matching edges. deba@2051: /// deba@2051: /// Gives back the number of the matching edges. deba@2051: int matchingSize() const { deba@2051: return matching_size; deba@2051: } deba@2051: deba@2051: /// @} deba@2051: deba@2051: private: deba@2051: deba@2051: void initStructures() { deba@2051: if (!_heap_cross_ref) { deba@2051: local_heap_cross_ref = true; deba@2051: _heap_cross_ref = Traits::createHeapCrossRef(*graph); deba@2051: } deba@2051: if (!_heap) { deba@2051: local_heap = true; deba@2051: _heap = Traits::createHeap(*_heap_cross_ref); deba@2051: } deba@2051: } deba@2051: deba@2051: void destroyStructures() { deba@2051: if (local_heap_cross_ref) delete _heap_cross_ref; deba@2051: if (local_heap) delete _heap; deba@2051: } deba@2051: deba@2051: deba@2051: private: deba@2051: deba@2051: const BpUGraph *graph; deba@2051: const CostMap* cost; deba@2051: deba@2051: ANodeMatchingMap anode_matching; deba@2051: BNodeMatchingMap bnode_matching; deba@2051: deba@2051: ANodePotentialMap anode_potential; deba@2051: BNodePotentialMap bnode_potential; deba@2051: deba@2051: Value matching_cost; deba@2051: int matching_size; deba@2051: deba@2051: HeapCrossRef *_heap_cross_ref; deba@2051: bool local_heap_cross_ref; deba@2051: deba@2051: Heap *_heap; deba@2051: bool local_heap; deba@2040: deba@2040: }; deba@2040: deba@2058: /// \ingroup matching deba@2058: /// deba@2058: /// \brief Minimum cost maximum cardinality bipartite matching deba@2058: /// deba@2058: /// This function calculates the minimum cost matching of the maximum deba@2058: /// cardinality matchings of a bipartite graph. It gives back the matching deba@2058: /// in an undirected edge map. deba@2058: /// deba@2058: /// \param graph The bipartite graph. deba@2058: /// \param cost The undirected edge map which contains the costs. deba@2058: /// \retval matching The undirected edge map which will be set to deba@2058: /// the matching. deba@2058: /// \return The cost of the matching. deba@2058: template deba@2058: typename CostMap::Value deba@2058: minCostMaxBipartiteMatching(const BpUGraph& graph, deba@2058: const CostMap& cost, deba@2058: MatchingMap& matching) { deba@2058: MinCostMaxBipartiteMatching deba@2058: bpmatching(graph, cost); deba@2058: bpmatching.run(); deba@2058: bpmatching.matching(matching); deba@2058: return bpmatching.matchingCost(); deba@2058: } deba@2058: deba@2040: } deba@2040: deba@2040: #endif