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