/* -*- C++ -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library

 *

 * Copyright (C) 2003-2008

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

#define LEMON_PR_BIPARTITE_MATCHING

#include <lemon/graph_utils.h>

#include <lemon/iterable_maps.h>

#include <iostream>

#include <queue>

#include <lemon/elevator.h>

///\ingroup matching

///\file

///\brief Push-prelabel maximum matching algorithms in bipartite graphs.

///

namespace lemon {

  ///Max cardinality matching algorithm based on push-relabel principle

  ///\ingroup matching

  ///Bipartite Max Cardinality Matching algorithm. This class uses the

  ///push-relabel principle which in several cases has better runtime

  ///performance than the augmenting path solutions.

  ///

  ///\author Alpar Juttner

  template<class Graph>

  class PrBipartiteMatching {

    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;

    const Graph &_g;

    int _node_num;

    int _matching_size;

    int _empty_level;

    typename Graph::template ANodeMap<typename Graph::UEdge> _matching;

    Elevator<Graph,typename Graph::BNode> _levels;

    typename Graph::template BNodeMap<int> _cov;

  public:

    /// Constructor

    /// Constructor

    ///

    PrBipartiteMatching(const Graph &g) :

      _g(g),

      _node_num(countBNodes(g)),

      _matching(g),

      _levels(g,_node_num),

      _cov(g,0)

    {

    }

    /// \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() and then one variant of the start()

    /// member. 

    /// @{

    ///Initialize the data structures

    ///This function constructs a prematching first, which is a

    ///regular matching on the A-side of the graph, but on the B-side

    ///each node could cover more matching edges. After that, the

    ///B-nodes which multiple matched, will be pushed into the lowest

    ///level of the Elevator. The remaning B-nodes will be pushed to

    ///the consequent levels respect to a Bfs on following graph: the

    ///nodes are the B-nodes of the original bipartite graph and two

    ///nodes are adjacent if a node can pass over a matching edge to

    ///an other node. The source of the Bfs are the lowest level

    ///nodes. Last, the reached B-nodes without covered matching edge

    ///becomes active.

    void init() {

      _matching_size=0;

      _empty_level=_node_num;

      for(ANodeIt n(_g);n!=INVALID;++n)

	if((_matching[n]=IncEdgeIt(_g,n))!=INVALID)

	  ++_cov[_g.bNode(_matching[n])];

      std::queue<Node> q;

      _levels.initStart();

      for(BNodeIt n(_g);n!=INVALID;++n)

	if(_cov[n]>1) {

	  _levels.initAddItem(n);

	  q.push(n);

	}

      int hlev=0;

      while(!q.empty()) {

	Node n=q.front();

	q.pop();

	int nlev=_levels[n]+1;

	for(IncEdgeIt e(_g,n);e!=INVALID;++e) {

	  Node m=_g.runningNode(e);

	  if(e==_matching[m]) {

	    for(IncEdgeIt f(_g,m);f!=INVALID;++f) {

	      Node r=_g.runningNode(f);

	      if(_levels[r]>nlev) {

		for(;nlev>hlev;hlev++)

		  _levels.initNewLevel();

		_levels.initAddItem(r);

		q.push(r);

	      }

	    }

	  }

	}

      }

      _levels.initFinish();

      for(BNodeIt n(_g);n!=INVALID;++n)

	if(_cov[n]<1&&_levels[n]<_node_num)

	  _levels.activate(n);

    }

    ///Start the main phase of the algorithm

    ///This algorithm calculates the maximum matching with the

    ///push-relabel principle. This function should be called just

    ///after the init() function which already set the initial

    ///prematching, the level function on the B-nodes and the active,

    ///ie. unmatched B-nodes.

    ///

    ///The algorithm always takes highest active B-node, and it try to

    ///find a B-node which is eligible to pass over one of it's

    ///matching edge. This condition holds when the B-node is one

    ///level lower, and the opposite node of it's matching edge is

    ///adjacent to the highest active node. In this case the current

    ///node steals the matching edge and becomes inactive. If there is

    ///not eligible node then the highest active node should be lift

    ///to the next proper level.

    ///

    ///The nodes should not lift higher than the number of the

    ///B-nodes, if a node reach this level it remains unmatched. If

    ///during the execution one level becomes empty the nodes above it

    ///can be deactivated and lift to the highest level.

    void start() {

      Node act;

      Node bact=INVALID;

      Node last_activated=INVALID;

      while((act=_levels.highestActive())!=INVALID) {

	last_activated=INVALID;

	int actlevel=_levels[act];

	UEdge bedge=INVALID;

	int nlevel=_node_num;

	{

	  int nnlevel;

	  for(IncEdgeIt tbedge(_g,act);

	      tbedge!=INVALID && nlevel>=actlevel;

	      ++tbedge)

	    if((nnlevel=_levels[_g.bNode(_matching[_g.runningNode(tbedge)])])<

	       nlevel)

	      {

		nlevel=nnlevel;

		bedge=tbedge;

	      }

	}

	if(nlevel<_node_num) {

	  if(nlevel>=actlevel)

	    _levels.liftHighestActive(nlevel+1);

	  bact=_g.bNode(_matching[_g.aNode(bedge)]);

	  if(--_cov[bact]<1) {

	    _levels.activate(bact);

	    last_activated=bact;

	  }

	  _matching[_g.aNode(bedge)]=bedge;

	  _cov[act]=1;

	  _levels.deactivate(act);

	}

	else {

	  _levels.liftHighestActiveToTop();

	}

	if(_levels.emptyLevel(actlevel))

	  _levels.liftToTop(actlevel);

      }

      for(ANodeIt n(_g);n!=INVALID;++n) {

	if (_matching[n]==INVALID)continue;

	if (_cov[_g.bNode(_matching[n])]>1) {

	  _cov[_g.bNode(_matching[n])]--;

	  _matching[n]=INVALID;

	} else {

	  ++_matching_size;

	}

      }

    }

    ///Start the algorithm to find a perfect matching

    ///This function is close to identical to the simple start()

    ///member function but it calculates just perfect matching.

    ///However, the perfect property is only checked on the B-side of

    ///the graph

    ///

    ///The main difference between the two function is the handling of

    ///the empty levels. The simple start() function let the nodes

    ///above the empty levels unmatched while this variant if it find

    ///an empty level immediately terminates and gives back false

    ///return value.

    bool startPerfect() {

      Node act;

      Node bact=INVALID;

      Node last_activated=INVALID;

      while((act=_levels.highestActive())!=INVALID) {

	last_activated=INVALID;

	int actlevel=_levels[act];

	UEdge bedge=INVALID;

	int nlevel=_node_num;

	{

	  int nnlevel;

	  for(IncEdgeIt tbedge(_g,act);

	      tbedge!=INVALID && nlevel>=actlevel;

	      ++tbedge)

	    if((nnlevel=_levels[_g.bNode(_matching[_g.runningNode(tbedge)])])<

	       nlevel)

	      {

		nlevel=nnlevel;

		bedge=tbedge;

	      }

	}

	if(nlevel<_node_num) {

	  if(nlevel>=actlevel)

	    _levels.liftHighestActive(nlevel+1);

	  bact=_g.bNode(_matching[_g.aNode(bedge)]);

	  if(--_cov[bact]<1) {

	    _levels.activate(bact);

	    last_activated=bact;

	  }

	  _matching[_g.aNode(bedge)]=bedge;

	  _cov[act]=1;

	  _levels.deactivate(act);

	}

	else {

	  _levels.liftHighestActiveToTop();

	}

	if(_levels.emptyLevel(actlevel))

	  _empty_level=actlevel;

	  return false;

      }

      _matching_size = _node_num;

      return true;

    }

    ///Runs the algorithm

    ///Just a shortcut for the next code:

    ///\code

    /// init();

    /// start();

    ///\endcode

    void run() {

      init();

      start();

    }

    ///Runs the algorithm to find a perfect matching

    ///Just a shortcut for the next code:

    ///\code

    /// init();

    /// startPerfect();

    ///\endcode

    ///

    ///\note If the two nodesets of the graph have different size then

    ///this algorithm checks the perfect property on the B-side.

    bool runPerfect() {

      init();

      return startPerfect();

    }

    ///Runs the algorithm to find a perfect matching

    ///Just a shortcut for the next code:

    ///\code

    /// init();

    /// startPerfect();

    ///\endcode

    ///

    ///\note It checks that the size of the two nodesets are equal.

    bool checkedRunPerfect() {

      if (countANodes(_g) != _node_num) return false;

      init();

      return startPerfect();

    }

    ///@}

    /// \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.

    ///@{

    ///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& mm) const {

      for (ANodeIt n(_g);n!=INVALID;++n) {

        if (_matching[n]!=INVALID) mm.set(_matching[n],true);

      }

      return _matching_size;

    }

    ///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<class MatchingMap>

    int matching(MatchingMap& mm) const {

      for (UEdgeIt e(_g);e!=INVALID;++e) {

        mm.set(e,e==_matching[_g.aNode(e)]);

      }

      return _matching_size;

    }

    ///Gives back the matching in an ANodeMap.

    ///Gives back the matching in an ANodeMap. The parameter should

    ///be a write ANodeMap of UEdge values.

    ///\return The number of the matching edges.

    template<class MatchingMap>

    int aMatching(MatchingMap& mm) const {

      for (ANodeIt n(_g);n!=INVALID;++n) {

        mm.set(n,_matching[n]);

      }

      return _matching_size;

    }

    ///Gives back the matching in a BNodeMap.

    ///Gives back the matching in a BNodeMap. The parameter should

    ///be a write BNodeMap of UEdge values.

    ///\return The number of the matching edges.

    template<class MatchingMap>

    int bMatching(MatchingMap& mm) const {

      for (BNodeIt n(_g);n!=INVALID;++n) {

        mm.set(n,INVALID);

      }

      for (ANodeIt n(_g);n!=INVALID;++n) {

        if (_matching[n]!=INVALID)

	  mm.set(_g.bNode(_matching[n]),_matching[n]);

      }

      return _matching_size;

    }

    ///Returns true if the given uedge is in the matching.

    ///It returns true if the given uedge is in the matching.

    ///

    bool matchingEdge(const UEdge& e) const {

      return _matching[_g.aNode(e)]==e;

    }

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

    ///\note If the parameter node is a B-node then the running time is

    ///propotional to the degree of the node.

    UEdge matchingEdge(const Node& n) const {

      if (_g.aNode(n)) {

        return _matching[n];

      } else {

	for (IncEdgeIt e(_g,n);e!=INVALID;++e)

	  if (e==_matching[_g.aNode(e)]) return e;

	return INVALID;

      }

    }

    ///Gives back the number of the matching edges.

    ///Gives back the number of the matching edges.

    int matchingSize() const {

      return _matching_size;

    }

    ///Gives back a barrier on the A-nodes

    ///The barrier is s subset of the nodes on the same side of the

    ///graph. If we tried to find a perfect matching and it failed

    ///then the barrier size will be greater than its neighbours. If

    ///the maximum matching searched then the barrier size minus its

    ///neighbours will be exactly the unmatched nodes on the A-side.

    ///\retval bar A WriteMap on the ANodes with bool value.

    template<class BarrierMap>

    void aBarrier(BarrierMap &bar) const 

    {

      for(ANodeIt n(_g);n!=INVALID;++n)

	bar.set(n,_matching[n]==INVALID ||

	  _levels[_g.bNode(_matching[n])]<_empty_level);  

    }  

    ///Gives back a barrier on the B-nodes

    ///The barrier is s subset of the nodes on the same side of the

    ///graph. If we tried to find a perfect matching and it failed

    ///then the barrier size will be greater than its neighbours. If

    ///the maximum matching searched then the barrier size minus its

    ///neighbours will be exactly the unmatched nodes on the B-side.

    ///\retval bar A WriteMap on the BNodes with bool value.

    template<class BarrierMap>

    void bBarrier(BarrierMap &bar) const

    {

      for(BNodeIt n(_g);n!=INVALID;++n) bar.set(n,_levels[n]>=_empty_level);  

    }

    ///Returns a minimum covering of the nodes.

    ///The minimum covering set problem is the dual solution of the

    ///maximum bipartite matching. It provides a solution for this

    ///problem what is proof of the optimality of the matching.

    ///\param covering NodeMap of bool values, the nodes of the cover

    ///set will set true while the others false.  

    ///\return The size of the cover set.

    ///\note This function can be called just after the algorithm have

    ///already found a matching. 

    template<class CoverMap>

    int coverSet(CoverMap& covering) const {

      int ret=0;

      for(BNodeIt n(_g);n!=INVALID;++n) {

	if (_levels[n]<_empty_level) { covering.set(n,true); ++ret; }

	else covering.set(n,false);

      }

      for(ANodeIt n(_g);n!=INVALID;++n) {

	if (_matching[n]!=INVALID &&

	    _levels[_g.bNode(_matching[n])]>=_empty_level) 

	  { covering.set(n,true); ++ret; }

	else covering.set(n,false);

      }

      return ret;

    }

    /// @}

  };

  ///Maximum cardinality of the matchings in a bipartite graph

  ///\ingroup matching

  ///This function finds the maximum cardinality of the matchings

  ///in a bipartite graph \c g.

  ///\param g An undirected bipartite graph.

  ///\return The cardinality of the maximum matching.

  ///

  ///\note The the implementation is based

  ///on the push-relabel principle.

  template<class Graph>

  int prBipartiteMatching(const Graph &g)

  {

    PrBipartiteMatching<Graph> bpm(g);

    return bpm.matchingSize();

  }

  ///Maximum cardinality matching in a bipartite graph

  ///\ingroup matching

  ///This function finds a maximum cardinality matching

  ///in a bipartite graph \c g.

  ///\param g An undirected bipartite graph.

  ///\retval matching A write ANodeMap of value type \c UEdge.

  /// The found edges will be returned in this map,

  /// i.e. for an \c ANode \c n the edge <tt>matching[n]</tt> is the one

  /// that covers the node \c n.

  ///\return The cardinality of the maximum matching.

  ///

  ///\note The the implementation is based

  ///on the push-relabel principle.

  template<class Graph,class MT>

  int prBipartiteMatching(const Graph &g,MT &matching) 

  {

    PrBipartiteMatching<Graph> bpm(g);

    bpm.run();

    bpm.aMatching(matching);

    return bpm.matchingSize();

  }

  ///Maximum cardinality matching in a bipartite graph

  ///\ingroup matching

  ///This function finds a maximum cardinality matching

  ///in a bipartite graph \c g.

  ///\param g An undirected bipartite graph.

  ///\retval matching A write ANodeMap of value type \c UEdge.

  /// The found edges will be returned in this map,

  /// i.e. for an \c ANode \c n the edge <tt>matching[n]</tt> is the one

  /// that covers the node \c n.

  ///\retval barrier A \c bool WriteMap on the BNodes. The map will be set

  /// exactly once for each BNode. The nodes with \c true value represent

  /// a barrier \e B, i.e. the cardinality of \e B minus the number of its

  /// neighbor is equal to the number of the <tt>BNode</tt>s minus the

  /// cardinality of the maximum matching.

  ///\return The cardinality of the maximum matching.

  ///

  ///\note The the implementation is based

  ///on the push-relabel principle.

  template<class Graph,class MT, class GT>

  int prBipartiteMatching(const Graph &g,MT &matching,GT &barrier) 

  {

    PrBipartiteMatching<Graph> bpm(g);

    bpm.run();

    bpm.aMatching(matching);

    bpm.bBarrier(barrier);

    return bpm.matchingSize();

  }  

  ///Perfect matching in a bipartite graph

  ///\ingroup matching

  ///This function checks whether the bipartite graph \c g

  ///has a perfect matching.

  ///\param g An undirected bipartite graph.

  ///\return \c true iff \c g has a perfect matching.

  ///

  ///\note The the implementation is based

  ///on the push-relabel principle.

  template<class Graph>

  bool prPerfectBipartiteMatching(const Graph &g)

  {

    PrBipartiteMatching<Graph> bpm(g);

    return bpm.runPerfect();

  }

  ///Perfect matching in a bipartite graph

  ///\ingroup matching

  ///This function finds a perfect matching in a bipartite graph \c g.

  ///\param g An undirected bipartite graph.

  ///\retval matching A write ANodeMap of value type \c UEdge.

  /// The found edges will be returned in this map,

  /// i.e. for an \c ANode \c n the edge <tt>matching[n]</tt> is the one

  /// that covers the node \c n.

  /// The values are unchanged if the graph

  /// has no perfect matching.

  ///\return \c true iff \c g has a perfect matching.

  ///

  ///\note The the implementation is based

  ///on the push-relabel principle.

  template<class Graph,class MT>

  bool prPerfectBipartiteMatching(const Graph &g,MT &matching) 

  {

    PrBipartiteMatching<Graph> bpm(g);

    bool ret = bpm.checkedRunPerfect();

    if (ret) bpm.aMatching(matching);

    return ret;

  }

  ///Perfect matching in a bipartite graph

  ///\ingroup matching

  ///This function finds a perfect matching in a bipartite graph \c g.

  ///\param g An undirected bipartite graph.

  ///\retval matching A write ANodeMap of value type \c UEdge.

  /// The found edges will be returned in this map,

  /// i.e. for an \c ANode \c n the edge <tt>matching[n]</tt> is the one

  /// that covers the node \c n.

  /// The values are unchanged if the graph

  /// has no perfect matching.

  ///\retval barrier A \c bool WriteMap on the BNodes. The map will only

  /// be set if \c g has no perfect matching. In this case it is set 

  /// exactly once for each BNode. The nodes with \c true value represent

  /// a barrier, i.e. a subset \e B a of BNodes with the property that

  /// the cardinality of \e B is greater than the number of its neighbors.

  ///\return \c true iff \c g has a perfect matching.

  ///

  ///\note The the implementation is based

  ///on the push-relabel principle.

  template<class Graph,class MT, class GT>

  bool prPerfectBipartiteMatching(const Graph &g,MT &matching,GT &barrier) 

  {

    PrBipartiteMatching<Graph> bpm(g);

    bool ret=bpm.checkedRunPerfect();

    if(ret)

      bpm.aMatching(matching);

    else

      bpm.bBarrier(barrier);

    return ret;

  }  

}

#endif