/* -*- 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 <lemon/bpugraph_adaptor.h>

 #include <lemon/bfs.h>

 #include <iostream>

 ///\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 <typename BpUGraph>

   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<UEdge> ANodeMatchingMap;

     typedef typename BpUGraph::template BNodeMap<UEdge> 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 <typename MatchingMap>

     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 <typename MatchingMap>

     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<bool> areached(*graph, false);

       typename Graph::template BNodeMap<bool> breached(*graph, false);

       typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);

       std::vector<Node> 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<Node> 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<bool> 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<bool> areached(*graph, false);

       typename Graph::template BNodeMap<bool> breached(*graph, false);

       typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);

       std::vector<Node> 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<Node> 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 <typename CoverMap>

     int coverSet(CoverMap& covering) {

       typename Graph::template ANodeMap<bool> areached(*graph, false);

       typename Graph::template BNodeMap<bool> breached(*graph, false);

       std::vector<Node> queue;

       for (ANodeIt it(*graph); it != INVALID; ++it) {

         if (anode_matching[it] == INVALID) {

           queue.push_back(it);

         }

       }

       while (!queue.empty()) {

         std::vector<Node> 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 <typename MatchingMap>

     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 <typename MatchingMap>

     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