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@2040: #include <lemon/bpugraph_adaptor.h>
deba@2040: #include <lemon/bfs.h>
deba@2040: 
deba@2040: #include <iostream>
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@2040:   /// the Hopcroft-Karp algorithm wich has \f$ O(e\sqrt{n}) \f$ time
deba@2040:   /// complexity.
deba@2040:   template <typename BpUGraph>
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<UEdge> ANodeMatchingMap;
deba@2040:     typedef typename BpUGraph::template BNodeMap<UEdge> 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@2040:       matching_value = 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@2040:       matching_value = 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@2040:             ++matching_value;
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 <typename MatchingMap>
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@2040:       matching_value = 0;
deba@2040:       for (UEdgeIt it(*graph); it != INVALID; ++it) {
deba@2040:         if (matching[it]) {
deba@2040:           ++matching_value;
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 <typename MatchingMap>
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@2040:       matching_value = 0;
deba@2040:       for (UEdgeIt it(*graph); it != INVALID; ++it) {
deba@2040:         if (matching[it]) {
deba@2040:           ++matching_value;
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<bool> areached(*graph, false);
deba@2040:       typename Graph::template BNodeMap<bool> breached(*graph, false);
deba@2040: 
deba@2040:       typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
deba@2040: 
deba@2040:       std::vector<Node> 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<Node> 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<bool> 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@2040:           ++matching_value;
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<bool> areached(*graph, false);
deba@2040:       typename Graph::template BNodeMap<bool> breached(*graph, false);
deba@2040: 
deba@2040:       typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
deba@2040: 
deba@2040:       std::vector<Node> 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<Node> 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@2040:               ++matching_value;
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@2040:       init();
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 <typename CoverMap>
deba@2040:     int coverSet(CoverMap& covering) {
deba@2040: 
deba@2040:       typename Graph::template ANodeMap<bool> areached(*graph, false);
deba@2040:       typename Graph::template BNodeMap<bool> breached(*graph, false);
deba@2040:       
deba@2040:       std::vector<Node> 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<Node> 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 <typename MatchingMap>
deba@2040:     int quickMatching(MatchingMap& matching) {
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@2040:       return matching_value;
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 <typename MatchingMap>
deba@2040:     int matching(MatchingMap& matching) {
deba@2040:       for (UEdgeIt it(*graph); it != INVALID; ++it) {
deba@2040:         matching[it] = it == anode_matching[graph->aNode(it)];
deba@2040:       }
deba@2040:       return matching_value;
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@2040:     bool matchingEdge(const UEdge& edge) {
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@2040:     UEdge matchingEdge(const Node& node) {
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@2040:     int matchingValue() const {
deba@2040:       return matching_value;
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@2040:     int matching_value;
deba@2040:   
deba@2040:   };
deba@2040: 
deba@2040: }
deba@2040: 
deba@2040: #endif