1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/bipartite_matching.h Fri Apr 07 09:51:23 2006 +0000
1.3 @@ -0,0 +1,466 @@
1.4 +/* -*- C++ -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library
1.7 + *
1.8 + * Copyright (C) 2003-2006
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#ifndef LEMON_BIPARTITE_MATCHING
1.23 +#define LEMON_BIPARTITE_MATCHING
1.24 +
1.25 +#include <lemon/bpugraph_adaptor.h>
1.26 +#include <lemon/bfs.h>
1.27 +
1.28 +#include <iostream>
1.29 +
1.30 +///\ingroup matching
1.31 +///\file
1.32 +///\brief Maximum matching algorithms in bipartite graphs.
1.33 +
1.34 +namespace lemon {
1.35 +
1.36 + /// \ingroup matching
1.37 + ///
1.38 + /// \brief Bipartite Max Cardinality Matching algorithm
1.39 + ///
1.40 + /// Bipartite Max Cardinality Matching algorithm. This class implements
1.41 + /// the Hopcroft-Karp algorithm wich has \f$ O(e\sqrt{n}) \f$ time
1.42 + /// complexity.
1.43 + template <typename BpUGraph>
1.44 + class MaxBipartiteMatching {
1.45 + protected:
1.46 +
1.47 + typedef BpUGraph Graph;
1.48 +
1.49 + typedef typename Graph::Node Node;
1.50 + typedef typename Graph::ANodeIt ANodeIt;
1.51 + typedef typename Graph::BNodeIt BNodeIt;
1.52 + typedef typename Graph::UEdge UEdge;
1.53 + typedef typename Graph::UEdgeIt UEdgeIt;
1.54 + typedef typename Graph::IncEdgeIt IncEdgeIt;
1.55 +
1.56 + typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
1.57 + typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
1.58 +
1.59 +
1.60 + public:
1.61 +
1.62 + /// \brief Constructor.
1.63 + ///
1.64 + /// Constructor of the algorithm.
1.65 + MaxBipartiteMatching(const BpUGraph& _graph)
1.66 + : anode_matching(_graph), bnode_matching(_graph), graph(&_graph) {}
1.67 +
1.68 + /// \name Execution control
1.69 + /// The simplest way to execute the algorithm is to use
1.70 + /// one of the member functions called \c run().
1.71 + /// \n
1.72 + /// If you need more control on the execution,
1.73 + /// first you must call \ref init() or one alternative for it.
1.74 + /// Finally \ref start() will perform the matching computation or
1.75 + /// with step-by-step execution you can augment the solution.
1.76 +
1.77 + /// @{
1.78 +
1.79 + /// \brief Initalize the data structures.
1.80 + ///
1.81 + /// It initalizes the data structures and creates an empty matching.
1.82 + void init() {
1.83 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.84 + anode_matching[it] = INVALID;
1.85 + }
1.86 + for (BNodeIt it(*graph); it != INVALID; ++it) {
1.87 + bnode_matching[it] = INVALID;
1.88 + }
1.89 + matching_value = 0;
1.90 + }
1.91 +
1.92 + /// \brief Initalize the data structures.
1.93 + ///
1.94 + /// It initalizes the data structures and creates a greedy
1.95 + /// matching. From this matching sometimes it is faster to get
1.96 + /// the matching than from the initial empty matching.
1.97 + void greedyInit() {
1.98 + matching_value = 0;
1.99 + for (BNodeIt it(*graph); it != INVALID; ++it) {
1.100 + bnode_matching[it] = INVALID;
1.101 + }
1.102 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.103 + anode_matching[it] = INVALID;
1.104 + for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
1.105 + if (bnode_matching[graph->bNode(jt)] == INVALID) {
1.106 + anode_matching[it] = jt;
1.107 + bnode_matching[graph->bNode(jt)] = jt;
1.108 + ++matching_value;
1.109 + break;
1.110 + }
1.111 + }
1.112 + }
1.113 + }
1.114 +
1.115 + /// \brief Initalize the data structures with an initial matching.
1.116 + ///
1.117 + /// It initalizes the data structures with an initial matching.
1.118 + template <typename MatchingMap>
1.119 + void matchingInit(const MatchingMap& matching) {
1.120 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.121 + anode_matching[it] = INVALID;
1.122 + }
1.123 + for (BNodeIt it(*graph); it != INVALID; ++it) {
1.124 + bnode_matching[it] = INVALID;
1.125 + }
1.126 + matching_value = 0;
1.127 + for (UEdgeIt it(*graph); it != INVALID; ++it) {
1.128 + if (matching[it]) {
1.129 + ++matching_value;
1.130 + anode_matching[graph->aNode(it)] = it;
1.131 + bnode_matching[graph->bNode(it)] = it;
1.132 + }
1.133 + }
1.134 + }
1.135 +
1.136 + /// \brief Initalize the data structures with an initial matching.
1.137 + ///
1.138 + /// It initalizes the data structures with an initial matching.
1.139 + /// \return %True when the given map contains really a matching.
1.140 + template <typename MatchingMap>
1.141 + void checkedMatchingInit(const MatchingMap& matching) {
1.142 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.143 + anode_matching[it] = INVALID;
1.144 + }
1.145 + for (BNodeIt it(*graph); it != INVALID; ++it) {
1.146 + bnode_matching[it] = INVALID;
1.147 + }
1.148 + matching_value = 0;
1.149 + for (UEdgeIt it(*graph); it != INVALID; ++it) {
1.150 + if (matching[it]) {
1.151 + ++matching_value;
1.152 + if (anode_matching[graph->aNode(it)] != INVALID) {
1.153 + return false;
1.154 + }
1.155 + anode_matching[graph->aNode(it)] = it;
1.156 + if (bnode_matching[graph->aNode(it)] != INVALID) {
1.157 + return false;
1.158 + }
1.159 + bnode_matching[graph->bNode(it)] = it;
1.160 + }
1.161 + }
1.162 + return false;
1.163 + }
1.164 +
1.165 + /// \brief An augmenting phase of the Hopcroft-Karp algorithm
1.166 + ///
1.167 + /// It runs an augmenting phase of the Hopcroft-Karp
1.168 + /// algorithm. The phase finds maximum count of edge disjoint
1.169 + /// augmenting paths and augments on these paths. The algorithm
1.170 + /// consists at most of \f$ O(\sqrt{n}) \f$ phase and one phase is
1.171 + /// \f$ O(e) \f$ long.
1.172 + bool augment() {
1.173 +
1.174 + typename Graph::template ANodeMap<bool> areached(*graph, false);
1.175 + typename Graph::template BNodeMap<bool> breached(*graph, false);
1.176 +
1.177 + typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
1.178 +
1.179 + std::vector<Node> queue, bqueue;
1.180 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.181 + if (anode_matching[it] == INVALID) {
1.182 + queue.push_back(it);
1.183 + areached[it] = true;
1.184 + }
1.185 + }
1.186 +
1.187 + bool success = false;
1.188 +
1.189 + while (!success && !queue.empty()) {
1.190 + std::vector<Node> newqueue;
1.191 + for (int i = 0; i < (int)queue.size(); ++i) {
1.192 + Node anode = queue[i];
1.193 + for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
1.194 + Node bnode = graph->bNode(jt);
1.195 + if (breached[bnode]) continue;
1.196 + breached[bnode] = true;
1.197 + bpred[bnode] = jt;
1.198 + if (bnode_matching[bnode] == INVALID) {
1.199 + bqueue.push_back(bnode);
1.200 + success = true;
1.201 + } else {
1.202 + Node newanode = graph->aNode(bnode_matching[bnode]);
1.203 + if (!areached[newanode]) {
1.204 + areached[newanode] = true;
1.205 + newqueue.push_back(newanode);
1.206 + }
1.207 + }
1.208 + }
1.209 + }
1.210 + queue.swap(newqueue);
1.211 + }
1.212 +
1.213 + if (success) {
1.214 +
1.215 + typename Graph::template ANodeMap<bool> aused(*graph, false);
1.216 +
1.217 + for (int i = 0; i < (int)bqueue.size(); ++i) {
1.218 + Node bnode = bqueue[i];
1.219 +
1.220 + bool used = false;
1.221 +
1.222 + while (bnode != INVALID) {
1.223 + UEdge uedge = bpred[bnode];
1.224 + Node anode = graph->aNode(uedge);
1.225 +
1.226 + if (aused[anode]) {
1.227 + used = true;
1.228 + break;
1.229 + }
1.230 +
1.231 + bnode = anode_matching[anode] != INVALID ?
1.232 + graph->bNode(anode_matching[anode]) : INVALID;
1.233 +
1.234 + }
1.235 +
1.236 + if (used) continue;
1.237 +
1.238 + bnode = bqueue[i];
1.239 + while (bnode != INVALID) {
1.240 + UEdge uedge = bpred[bnode];
1.241 + Node anode = graph->aNode(uedge);
1.242 +
1.243 + bnode_matching[bnode] = uedge;
1.244 +
1.245 + bnode = anode_matching[anode] != INVALID ?
1.246 + graph->bNode(anode_matching[anode]) : INVALID;
1.247 +
1.248 + anode_matching[anode] = uedge;
1.249 +
1.250 + aused[anode] = true;
1.251 + }
1.252 + ++matching_value;
1.253 +
1.254 + }
1.255 + }
1.256 + return success;
1.257 + }
1.258 +
1.259 + /// \brief An augmenting phase of the Ford-Fulkerson algorithm
1.260 + ///
1.261 + /// It runs an augmenting phase of the Ford-Fulkerson
1.262 + /// algorithm. The phase finds only one augmenting path and
1.263 + /// augments only on this paths. The algorithm consists at most
1.264 + /// of \f$ O(n) \f$ simple phase and one phase is at most
1.265 + /// \f$ O(e) \f$ long.
1.266 + bool simpleAugment() {
1.267 +
1.268 + typename Graph::template ANodeMap<bool> areached(*graph, false);
1.269 + typename Graph::template BNodeMap<bool> breached(*graph, false);
1.270 +
1.271 + typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
1.272 +
1.273 + std::vector<Node> queue;
1.274 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.275 + if (anode_matching[it] == INVALID) {
1.276 + queue.push_back(it);
1.277 + areached[it] = true;
1.278 + }
1.279 + }
1.280 +
1.281 + while (!queue.empty()) {
1.282 + std::vector<Node> newqueue;
1.283 + for (int i = 0; i < (int)queue.size(); ++i) {
1.284 + Node anode = queue[i];
1.285 + for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
1.286 + Node bnode = graph->bNode(jt);
1.287 + if (breached[bnode]) continue;
1.288 + breached[bnode] = true;
1.289 + bpred[bnode] = jt;
1.290 + if (bnode_matching[bnode] == INVALID) {
1.291 + while (bnode != INVALID) {
1.292 + UEdge uedge = bpred[bnode];
1.293 + anode = graph->aNode(uedge);
1.294 +
1.295 + bnode_matching[bnode] = uedge;
1.296 +
1.297 + bnode = anode_matching[anode] != INVALID ?
1.298 + graph->bNode(anode_matching[anode]) : INVALID;
1.299 +
1.300 + anode_matching[anode] = uedge;
1.301 +
1.302 + }
1.303 + ++matching_value;
1.304 + return true;
1.305 + } else {
1.306 + Node newanode = graph->aNode(bnode_matching[bnode]);
1.307 + if (!areached[newanode]) {
1.308 + areached[newanode] = true;
1.309 + newqueue.push_back(newanode);
1.310 + }
1.311 + }
1.312 + }
1.313 + }
1.314 + queue.swap(newqueue);
1.315 + }
1.316 +
1.317 + return false;
1.318 + }
1.319 +
1.320 + /// \brief Starts the algorithm.
1.321 + ///
1.322 + /// Starts the algorithm. It runs augmenting phases until the optimal
1.323 + /// solution reached.
1.324 + void start() {
1.325 + while (augment()) {}
1.326 + }
1.327 +
1.328 + /// \brief Runs the algorithm.
1.329 + ///
1.330 + /// It just initalize the algorithm and then start it.
1.331 + void run() {
1.332 + init();
1.333 + start();
1.334 + }
1.335 +
1.336 + /// @}
1.337 +
1.338 + /// \name Query Functions
1.339 + /// The result of the %Matching algorithm can be obtained using these
1.340 + /// functions.\n
1.341 + /// Before the use of these functions,
1.342 + /// either run() or start() must be called.
1.343 +
1.344 + ///@{
1.345 +
1.346 + /// \brief Returns an minimum covering of the nodes.
1.347 + ///
1.348 + /// The minimum covering set problem is the dual solution of the
1.349 + /// maximum bipartite matching. It provides an solution for this
1.350 + /// problem what is proof of the optimality of the matching.
1.351 + /// \return The size of the cover set.
1.352 + template <typename CoverMap>
1.353 + int coverSet(CoverMap& covering) {
1.354 +
1.355 + typename Graph::template ANodeMap<bool> areached(*graph, false);
1.356 + typename Graph::template BNodeMap<bool> breached(*graph, false);
1.357 +
1.358 + std::vector<Node> queue;
1.359 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.360 + if (anode_matching[it] == INVALID) {
1.361 + queue.push_back(it);
1.362 + }
1.363 + }
1.364 +
1.365 + while (!queue.empty()) {
1.366 + std::vector<Node> newqueue;
1.367 + for (int i = 0; i < (int)queue.size(); ++i) {
1.368 + Node anode = queue[i];
1.369 + for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
1.370 + Node bnode = graph->bNode(jt);
1.371 + if (breached[bnode]) continue;
1.372 + breached[bnode] = true;
1.373 + if (bnode_matching[bnode] != INVALID) {
1.374 + Node newanode = graph->aNode(bnode_matching[bnode]);
1.375 + if (!areached[newanode]) {
1.376 + areached[newanode] = true;
1.377 + newqueue.push_back(newanode);
1.378 + }
1.379 + }
1.380 + }
1.381 + }
1.382 + queue.swap(newqueue);
1.383 + }
1.384 +
1.385 + int size = 0;
1.386 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.387 + covering[it] = !areached[it] && anode_matching[it] != INVALID;
1.388 + if (!areached[it] && anode_matching[it] != INVALID) {
1.389 + ++size;
1.390 + }
1.391 + }
1.392 + for (BNodeIt it(*graph); it != INVALID; ++it) {
1.393 + covering[it] = breached[it];
1.394 + if (breached[it]) {
1.395 + ++size;
1.396 + }
1.397 + }
1.398 + return size;
1.399 + }
1.400 +
1.401 + /// \brief Set true all matching uedge in the map.
1.402 + ///
1.403 + /// Set true all matching uedge in the map. It does not change the
1.404 + /// value mapped to the other uedges.
1.405 + /// \return The number of the matching edges.
1.406 + template <typename MatchingMap>
1.407 + int quickMatching(MatchingMap& matching) {
1.408 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.409 + if (anode_matching[it] != INVALID) {
1.410 + matching[anode_matching[it]] = true;
1.411 + }
1.412 + }
1.413 + return matching_value;
1.414 + }
1.415 +
1.416 + /// \brief Set true all matching uedge in the map and the others to false.
1.417 + ///
1.418 + /// Set true all matching uedge in the map and the others to false.
1.419 + /// \return The number of the matching edges.
1.420 + template <typename MatchingMap>
1.421 + int matching(MatchingMap& matching) {
1.422 + for (UEdgeIt it(*graph); it != INVALID; ++it) {
1.423 + matching[it] = it == anode_matching[graph->aNode(it)];
1.424 + }
1.425 + return matching_value;
1.426 + }
1.427 +
1.428 +
1.429 + /// \brief Return true if the given uedge is in the matching.
1.430 + ///
1.431 + /// It returns true if the given uedge is in the matching.
1.432 + bool matchingEdge(const UEdge& edge) {
1.433 + return anode_matching[graph->aNode(edge)] == edge;
1.434 + }
1.435 +
1.436 + /// \brief Returns the matching edge from the node.
1.437 + ///
1.438 + /// Returns the matching edge from the node. If there is not such
1.439 + /// edge it gives back \c INVALID.
1.440 + UEdge matchingEdge(const Node& node) {
1.441 + if (graph->aNode(node)) {
1.442 + return anode_matching[node];
1.443 + } else {
1.444 + return bnode_matching[node];
1.445 + }
1.446 + }
1.447 +
1.448 + /// \brief Gives back the number of the matching edges.
1.449 + ///
1.450 + /// Gives back the number of the matching edges.
1.451 + int matchingValue() const {
1.452 + return matching_value;
1.453 + }
1.454 +
1.455 + /// @}
1.456 +
1.457 + private:
1.458 +
1.459 + ANodeMatchingMap anode_matching;
1.460 + BNodeMatchingMap bnode_matching;
1.461 + const Graph *graph;
1.462 +
1.463 + int matching_value;
1.464 +
1.465 + };
1.466 +
1.467 +}
1.468 +
1.469 +#endif