# HG changeset patch # User deba # Date 1144403483 0 # Node ID c7bd55c0d8202776b2ff5010ee339a3cca82c7ba # Parent dacc4ce9474dddca29c8f1f8be55e3a7510121f1 Bipartite Graph Max Cardinality Matching (Hopcroft-Karp) Test for it Some BpUgraph improvments diff -r dacc4ce9474d -r c7bd55c0d820 lemon/Makefile.am --- a/lemon/Makefile.am Thu Apr 06 09:33:29 2006 +0000 +++ b/lemon/Makefile.am Fri Apr 07 09:51:23 2006 +0000 @@ -31,6 +31,8 @@ dfs.h \ bin_heap.h \ bpugraph_adaptor.h \ + bipartite_matching.h \ + bucket_heap.h \ color.h \ config.h \ counter.h \ @@ -53,7 +55,6 @@ iterable_maps.h \ johnson.h \ kruskal.h \ - bucket_heap.h \ list_graph.h \ lp.h \ lp_base.h \ diff -r dacc4ce9474d -r c7bd55c0d820 lemon/bipartite_matching.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lemon/bipartite_matching.h Fri Apr 07 09:51:23 2006 +0000 @@ -0,0 +1,466 @@ +/* -*- C++ -*- + * + * This file is a part of LEMON, a generic C++ optimization library + * + * Copyright (C) 2003-2006 + * 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 + +///\ingroup matching +///\file +///\brief Maximum matching algorithms in bipartite graphs. + +namespace lemon { + + /// \ingroup matching + /// + /// \brief Bipartite Max Cardinality Matching algorithm + /// + /// Bipartite Max Cardinality Matching algorithm. This class implements + /// the Hopcroft-Karp algorithm wich has \f$ O(e\sqrt{n}) \f$ time + /// complexity. + 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) + : anode_matching(_graph), bnode_matching(_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) { + anode_matching[it] = INVALID; + } + for (BNodeIt it(*graph); it != INVALID; ++it) { + bnode_matching[it] = INVALID; + } + matching_value = 0; + } + + /// \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() { + matching_value = 0; + for (BNodeIt it(*graph); it != INVALID; ++it) { + bnode_matching[it] = INVALID; + } + for (ANodeIt it(*graph); it != INVALID; ++it) { + anode_matching[it] = INVALID; + for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) { + if (bnode_matching[graph->bNode(jt)] == INVALID) { + anode_matching[it] = jt; + bnode_matching[graph->bNode(jt)] = jt; + ++matching_value; + break; + } + } + } + } + + /// \brief Initalize the data structures with an initial matching. + /// + /// It initalizes the data structures with an initial matching. + template + void matchingInit(const MatchingMap& matching) { + for (ANodeIt it(*graph); it != INVALID; ++it) { + anode_matching[it] = INVALID; + } + for (BNodeIt it(*graph); it != INVALID; ++it) { + bnode_matching[it] = INVALID; + } + matching_value = 0; + for (UEdgeIt it(*graph); it != INVALID; ++it) { + if (matching[it]) { + ++matching_value; + anode_matching[graph->aNode(it)] = it; + bnode_matching[graph->bNode(it)] = it; + } + } + } + + /// \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 + void checkedMatchingInit(const MatchingMap& matching) { + for (ANodeIt it(*graph); it != INVALID; ++it) { + anode_matching[it] = INVALID; + } + for (BNodeIt it(*graph); it != INVALID; ++it) { + bnode_matching[it] = INVALID; + } + matching_value = 0; + for (UEdgeIt it(*graph); it != INVALID; ++it) { + if (matching[it]) { + ++matching_value; + if (anode_matching[graph->aNode(it)] != INVALID) { + return false; + } + anode_matching[graph->aNode(it)] = it; + if (bnode_matching[graph->aNode(it)] != INVALID) { + return false; + } + bnode_matching[graph->bNode(it)] = it; + } + } + return false; + } + + /// \brief An augmenting phase of the Hopcroft-Karp algorithm + /// + /// It runs an augmenting phase of the Hopcroft-Karp + /// algorithm. The phase finds maximum count of 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() { + + typename Graph::template ANodeMap areached(*graph, false); + typename Graph::template BNodeMap breached(*graph, false); + + typename Graph::template BNodeMap bpred(*graph, INVALID); + + std::vector queue, bqueue; + for (ANodeIt it(*graph); it != INVALID; ++it) { + if (anode_matching[it] == INVALID) { + queue.push_back(it); + areached[it] = true; + } + } + + bool success = false; + + while (!success && !queue.empty()) { + std::vector newqueue; + for (int i = 0; i < (int)queue.size(); ++i) { + Node anode = queue[i]; + for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) { + Node bnode = graph->bNode(jt); + if (breached[bnode]) continue; + breached[bnode] = true; + bpred[bnode] = jt; + if (bnode_matching[bnode] == INVALID) { + bqueue.push_back(bnode); + success = true; + } else { + Node newanode = graph->aNode(bnode_matching[bnode]); + if (!areached[newanode]) { + areached[newanode] = true; + newqueue.push_back(newanode); + } + } + } + } + queue.swap(newqueue); + } + + if (success) { + + typename Graph::template ANodeMap aused(*graph, false); + + for (int i = 0; i < (int)bqueue.size(); ++i) { + Node bnode = bqueue[i]; + + bool used = false; + + while (bnode != INVALID) { + UEdge uedge = bpred[bnode]; + Node anode = graph->aNode(uedge); + + if (aused[anode]) { + used = true; + break; + } + + bnode = anode_matching[anode] != INVALID ? + graph->bNode(anode_matching[anode]) : INVALID; + + } + + if (used) continue; + + bnode = bqueue[i]; + while (bnode != INVALID) { + UEdge uedge = bpred[bnode]; + Node anode = graph->aNode(uedge); + + bnode_matching[bnode] = uedge; + + bnode = anode_matching[anode] != INVALID ? + graph->bNode(anode_matching[anode]) : INVALID; + + anode_matching[anode] = uedge; + + aused[anode] = true; + } + ++matching_value; + + } + } + return success; + } + + /// \brief An augmenting phase of the Ford-Fulkerson algorithm + /// + /// It runs an augmenting phase of the Ford-Fulkerson + /// algorithm. The phase finds only one augmenting path and + /// augments only on this paths. The algorithm consists at most + /// of \f$ O(n) \f$ simple phase and one phase is at most + /// \f$ O(e) \f$ long. + bool simpleAugment() { + + typename Graph::template ANodeMap areached(*graph, false); + typename Graph::template BNodeMap breached(*graph, false); + + typename Graph::template BNodeMap bpred(*graph, INVALID); + + std::vector queue; + for (ANodeIt it(*graph); it != INVALID; ++it) { + if (anode_matching[it] == INVALID) { + queue.push_back(it); + areached[it] = true; + } + } + + while (!queue.empty()) { + std::vector newqueue; + for (int i = 0; i < (int)queue.size(); ++i) { + Node anode = queue[i]; + for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) { + Node bnode = graph->bNode(jt); + if (breached[bnode]) continue; + breached[bnode] = true; + bpred[bnode] = jt; + if (bnode_matching[bnode] == INVALID) { + while (bnode != INVALID) { + UEdge uedge = bpred[bnode]; + anode = graph->aNode(uedge); + + bnode_matching[bnode] = uedge; + + bnode = anode_matching[anode] != INVALID ? + graph->bNode(anode_matching[anode]) : INVALID; + + anode_matching[anode] = uedge; + + } + ++matching_value; + return true; + } else { + Node newanode = graph->aNode(bnode_matching[bnode]); + if (!areached[newanode]) { + areached[newanode] = true; + newqueue.push_back(newanode); + } + } + } + } + queue.swap(newqueue); + } + + return false; + } + + /// \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 Returns an minimum covering of the nodes. + /// + /// The minimum covering set problem is the dual solution of the + /// maximum bipartite matching. It provides an 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) { + + typename Graph::template ANodeMap areached(*graph, false); + typename Graph::template BNodeMap breached(*graph, false); + + std::vector queue; + for (ANodeIt it(*graph); it != INVALID; ++it) { + if (anode_matching[it] == INVALID) { + queue.push_back(it); + } + } + + while (!queue.empty()) { + std::vector newqueue; + for (int i = 0; i < (int)queue.size(); ++i) { + Node anode = queue[i]; + for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) { + Node bnode = graph->bNode(jt); + if (breached[bnode]) continue; + breached[bnode] = true; + if (bnode_matching[bnode] != INVALID) { + Node newanode = graph->aNode(bnode_matching[bnode]); + if (!areached[newanode]) { + areached[newanode] = true; + newqueue.push_back(newanode); + } + } + } + } + queue.swap(newqueue); + } + + int size = 0; + for (ANodeIt it(*graph); it != INVALID; ++it) { + covering[it] = !areached[it] && anode_matching[it] != INVALID; + if (!areached[it] && anode_matching[it] != INVALID) { + ++size; + } + } + for (BNodeIt it(*graph); it != INVALID; ++it) { + covering[it] = breached[it]; + if (breached[it]) { + ++size; + } + } + return size; + } + + /// \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& matching) { + for (ANodeIt it(*graph); it != INVALID; ++it) { + if (anode_matching[it] != INVALID) { + matching[anode_matching[it]] = true; + } + } + return matching_value; + } + + /// \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& matching) { + for (UEdgeIt it(*graph); it != INVALID; ++it) { + matching[it] = it == anode_matching[graph->aNode(it)]; + } + return matching_value; + } + + + /// \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) { + 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) { + if (graph->aNode(node)) { + return anode_matching[node]; + } else { + return bnode_matching[node]; + } + } + + /// \brief Gives back the number of the matching edges. + /// + /// Gives back the number of the matching edges. + int matchingValue() const { + return matching_value; + } + + /// @} + + private: + + ANodeMatchingMap anode_matching; + BNodeMatchingMap bnode_matching; + const Graph *graph; + + int matching_value; + + }; + +} + +#endif diff -r dacc4ce9474d -r c7bd55c0d820 lemon/bpugraph_adaptor.h --- a/lemon/bpugraph_adaptor.h Thu Apr 06 09:33:29 2006 +0000 +++ b/lemon/bpugraph_adaptor.h Fri Apr 07 09:51:23 2006 +0000 @@ -103,6 +103,12 @@ Node source(const Edge& e) const { return graph->source(e); } Node target(const Edge& e) const { return graph->target(e); } + Node aNode(const UEdge& e) const { return graph->aNode(e); } + Node bNode(const UEdge& e) const { return graph->bNode(e); } + + bool aNode(const Node& n) const { return graph->aNode(n); } + bool bNode(const Node& n) const { return graph->bNode(n); } + typedef NodeNumTagIndicator NodeNumTag; int nodeNum() const { return graph->nodeNum(); } int aNodeNum() const { return graph->aNodeNum(); } @@ -122,7 +128,8 @@ return graph->findUEdge(source, target, prev); } - Node addNode() const { return graph->addNode(); } + Node addANode() const { return graph->addANode(); } + Node addBNode() const { return graph->addBNode(); } UEdge addEdge(const Node& source, const Node& target) const { return graph->addEdge(source, target); } @@ -281,6 +288,11 @@ }; /// \ingroup graph_adaptors + /// + /// \brief Trivial Bipartite Undirected Graph Adaptor + /// + /// Trivial Bipartite Undirected Graph Adaptor. It does not change + /// the functionality of the adapted graph. template class BpUGraphAdaptor : public BpUGraphAdaptorExtender< BpUGraphAdaptorBase<_BpUGraph> > { @@ -310,6 +322,17 @@ typedef typename Parent::Node Node; typedef typename Parent::BNode ANode; typedef typename Parent::ANode BNode; + typedef typename Parent::Edge Edge; + typedef typename Parent::UEdge UEdge; + + bool direction(const Edge& e) const { return !Parent::direction(e); } + Edge direct(const UEdge& e, bool b) const { return Parent::direct(e, !b); } + + Node aNode(const UEdge& e) const { return Parent::bNode(e); } + Node bNode(const UEdge& e) const { return Parent::aNode(e); } + + bool aNode(const Node& n) const { return Parent::bNode(n); } + bool bNode(const Node& n) const { return Parent::aNode(n); } void firstANode(Node& i) const { Parent::firstBNode(i); } void firstBNode(Node& i) const { Parent::firstANode(i); } @@ -407,6 +430,126 @@ }; + template + class MatchingBpUGraphAdaptorBase + : public BpUGraphAdaptorBase + { + public: + + typedef _BpUGraph Graph; + typedef _ANMatchingMap ANodeMatchingMap; + typedef _BNMatchingMap BNodeMatchingMap; + + typedef BpUGraphAdaptorBase Parent; + + protected: + + MatchingBpUGraphAdaptorBase() {} + + void setANodeMatchingMap(ANodeMatchingMap& _anode_matching) { + anode_matching = &_anode_matching; + } + + void setBNodeMatchingMap(BNodeMatchingMap& _bnode_matching) { + bnode_matching = &_bnode_matching; + } + + public: + + typedef typename Parent::Node Node; + typedef typename Parent::Edge Edge; + + void firstOut(Edge& edge, const Node& node) const { + if (Parent::aNode(node)) { + Parent::firstOut(edge, node); + if (edge != INVALID && (*anode_matching)[node] == edge) { + Parent::nextOut(edge); + } + } else { + edge = (*bnode_matching)[node] != INVALID ? + Parent::direct((*bnode_matching)[node], false) : INVALID; + } + } + + void nextOut(Edge& edge) const { + if (Parent::aNode(Parent::source(edge))) { + Parent::nextOut(edge); + if (edge != INVALID && (*anode_matching)[Parent::aNode(edge)] == edge) { + Parent::nextOut(edge); + } + } else { + edge = INVALID; + } + } + + void firstIn(Edge& edge, const Node& node) const { + if (Parent::aNode(node)) { + edge = (*bnode_matching)[node] != INVALID ? + Parent::direct((*anode_matching)[node], false) : INVALID; + } else { + Parent::firstIn(edge, node); + if (edge != INVALID && (*bnode_matching)[node] == edge) { + Parent::nextIn(edge); + } + } + } + + void nextIn(Edge& edge) const { + if (Parent::aNode(Parent::target(edge))) { + edge = INVALID; + } else { + Parent::nextIn(edge); + if (edge != INVALID && (*bnode_matching)[Parent::bNode(edge)] == edge) { + Parent::nextIn(edge); + } + } + } + + private: + ANodeMatchingMap* anode_matching; + BNodeMatchingMap* bnode_matching; + }; + + + /// \ingroup graph_adaptors + /// + /// \brief Bipartite graph adaptor to implement matching algorithms. + /// + /// Bipartite graph adaptor to implement matchings. It implements + /// the residual graph of the matching. + template + class MatchingBpUGraphAdaptor + : public BpUGraphAdaptorExtender< + MatchingBpUGraphAdaptorBase<_BpUGraph, _ANMatchingMap, _BNMatchingMap> > + { + public: + + typedef _BpUGraph BpUGraph; + typedef _BpUGraph Graph; + typedef _ANMatchingMap ANodeMatchingMap; + typedef _BNMatchingMap BNodeMatchingMap; + typedef BpUGraphAdaptorExtender< + MatchingBpUGraphAdaptorBase > + Parent; + + protected: + MatchingBpUGraphAdaptor() : Parent() {} + + public: + + MatchingBpUGraphAdaptor(const Graph& _graph, + ANodeMatchingMap& _anode_matching, + BNodeMatchingMap& _bnode_matching) { + setGraph(_graph); + setANodeMatchingMap(_anode_matching); + setBNodeMatchingMap(_bnode_matching); + } + + }; + } #endif diff -r dacc4ce9474d -r c7bd55c0d820 test/Makefile.am --- a/test/Makefile.am Thu Apr 06 09:33:29 2006 +0000 +++ b/test/Makefile.am Fri Apr 07 09:51:23 2006 +0000 @@ -15,6 +15,7 @@ all_pairs_shortest_path_test \ arborescence_test \ bfs_test \ + bipartite_matching_test \ counter_test \ dfs_test \ dijkstra_test \ @@ -55,6 +56,7 @@ all_pairs_shortest_path_test_SOURCES = all_pairs_shortest_path_test.cc arborescence_test_SOURCES = arborescence_test.cc bfs_test_SOURCES = bfs_test.cc +bipartite_matching_test_SOURCES = bipartite_matching_test.cc counter_test_SOURCES = counter_test.cc dfs_test_SOURCES = dfs_test.cc dijkstra_test_SOURCES = dijkstra_test.cc diff -r dacc4ce9474d -r c7bd55c0d820 test/bipartite_matching_test.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/bipartite_matching_test.cc Fri Apr 07 09:51:23 2006 +0000 @@ -0,0 +1,138 @@ +#include +#include +#include + +#include + +#include +#include + +#include + +#include "test_tools.h" + +using namespace std; +using namespace lemon; + +typedef ListBpUGraph Graph; +BPUGRAPH_TYPEDEFS(Graph); + +int main(int argc, const char *argv[]) { + Graph graph; + vector aNodes; + vector bNodes; + int n = argc > 1 ? atoi(argv[1]) : 100; + int m = argc > 2 ? atoi(argv[2]) : 100; + int e = argc > 3 ? atoi(argv[3]) : (int)((n+m) * log(n+m)); + + for (int i = 0; i < n; ++i) { + Node node = graph.addANode(); + aNodes.push_back(node); + } + for (int i = 0; i < m; ++i) { + Node node = graph.addBNode(); + bNodes.push_back(node); + } + for (int i = 0; i < e; ++i) { + Node aNode = aNodes[urandom(n)]; + Node bNode = bNodes[urandom(m)]; + graph.addEdge(aNode, bNode); + } + + { + MaxBipartiteMatching bpmatch(graph); + + bpmatch.run(); + + Graph::UEdgeMap mm(graph); + Graph::NodeMap cs(graph); + + check(bpmatch.coverSet(cs) == bpmatch.matching(mm), "INVALID PRIMAL-DUAL"); + + for (UEdgeIt it(graph); it != INVALID; ++it) { + check(cs[graph.aNode(it)] || cs[graph.bNode(it)], "INVALID DUAL"); + } + + for (ANodeIt it(graph); it != INVALID; ++it) { + int num = 0; + for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) { + if (mm[jt]) ++num; + } + check(num <= 1, "INVALID PRIMAL"); + } + } + + { + MaxBipartiteMatching bpmatch(graph); + + bpmatch.greedyInit(); + bpmatch.start(); + + Graph::UEdgeMap mm(graph); + Graph::NodeMap cs(graph); + + check(bpmatch.coverSet(cs) == bpmatch.matching(mm), "INVALID PRIMAL-DUAL"); + + for (UEdgeIt it(graph); it != INVALID; ++it) { + check(cs[graph.aNode(it)] || cs[graph.bNode(it)], "INVALID DUAL"); + } + + for (ANodeIt it(graph); it != INVALID; ++it) { + int num = 0; + for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) { + if (mm[jt]) ++num; + } + check(num <= 1, "INVALID PRIMAL"); + } + } + + { + MaxBipartiteMatching bpmatch(graph); + + bpmatch.greedyInit(); + while (bpmatch.simpleAugment()); + + Graph::UEdgeMap mm(graph); + Graph::NodeMap cs(graph); + + check(bpmatch.coverSet(cs) == bpmatch.matching(mm), "INVALID PRIMAL-DUAL"); + + for (UEdgeIt it(graph); it != INVALID; ++it) { + check(cs[graph.aNode(it)] || cs[graph.bNode(it)], "INVALID DUAL"); + } + + for (ANodeIt it(graph); it != INVALID; ++it) { + int num = 0; + for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) { + if (mm[jt]) ++num; + } + check(num <= 1, "INVALID PRIMAL"); + } + } + + { + SwapBpUGraphAdaptor swapgraph(graph); + MaxBipartiteMatching > bpmatch(swapgraph); + + bpmatch.run(); + + Graph::UEdgeMap mm(graph); + Graph::NodeMap cs(graph); + + check(bpmatch.coverSet(cs) == bpmatch.matching(mm), "INVALID PRIMAL-DUAL"); + + for (UEdgeIt it(graph); it != INVALID; ++it) { + check(cs[graph.aNode(it)] || cs[graph.bNode(it)], "INVALID DUAL"); + } + + for (ANodeIt it(graph); it != INVALID; ++it) { + int num = 0; + for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) { + if (mm[jt]) ++num; + } + check(num <= 1, "INVALID PRIMAL"); + } + } + + return 0; +}