/* -*- 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