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

#define LEMON_MAX_MATCHING_H

#include <queue>

#include <lemon/bits/invalid.h>

#include <lemon/unionfind.h>

#include <lemon/graph_utils.h>

///\ingroup matching

///\file

///\brief Maximum matching algorithm in undirected graph.

namespace lemon {

  /// \ingroup matching

  ///Edmonds' alternating forest maximum matching algorithm.

  ///This class provides Edmonds' alternating forest matching

  ///algorithm. The starting matching (if any) can be passed to the

  ///algorithm using read-in functions \ref readNMapNode, \ref

  ///readNMapEdge or \ref readEMapBool depending on the container. The

  ///resulting maximum matching can be attained by write-out functions

  ///\ref writeNMapNode, \ref writeNMapEdge or \ref writeEMapBool

  ///depending on the preferred container. 

  ///

  ///The dual side of a matching is a map of the nodes to

  ///MaxMatching::pos_enum, having values D, A and C showing the

  ///Gallai-Edmonds decomposition of the graph. The nodes in D induce

  ///a graph with factor-critical components, the nodes in A form the

  ///barrier, and the nodes in C induce a graph having a perfect

  ///matching. This decomposition can be attained by calling \ref

  ///writePos after running the algorithm. 

  ///

  ///\param Graph The undirected graph type the algorithm runs on.

  ///

  ///\author Jacint Szabo  

  template <typename Graph>

  class MaxMatching {

  protected:

    typedef typename Graph::Node Node;

    typedef typename Graph::Edge Edge;

    typedef typename Graph::UEdge UEdge;

    typedef typename Graph::UEdgeIt UEdgeIt;

    typedef typename Graph::NodeIt NodeIt;

    typedef typename Graph::IncEdgeIt IncEdgeIt;

    typedef typename Graph::template NodeMap<int> UFECrossRef;

    typedef UnionFindEnum<Node, UFECrossRef> UFE;

  public:

    ///Indicates the Gallai-Edmonds decomposition of the graph.

    ///Indicates the Gallai-Edmonds decomposition of the graph, which

    ///shows an upper bound on the size of a maximum matching. The

    ///nodes with pos_enum \c D induce a graph with factor-critical

    ///components, the nodes in \c A form the canonical barrier, and the

    ///nodes in \c C induce a graph having a perfect matching. 

    enum pos_enum {

      D=0,

      A=1,

      C=2

    }; 

  protected:

    static const int HEUR_density=2;

    const Graph& g;

    typename Graph::template NodeMap<Node> _mate;

    typename Graph::template NodeMap<pos_enum> position;

  public:

    MaxMatching(const Graph& _g) : g(_g), _mate(_g,INVALID), position(_g) {}

    ///Runs Edmonds' algorithm.

    ///Runs Edmonds' algorithm for sparse graphs (number of edges <

    ///2*number of nodes), and a heuristical Edmonds' algorithm with a

    ///heuristic of postponing shrinks for dense graphs. 

    void run() {

      if ( countUEdges(g) < HEUR_density*countNodes(g) ) {

	greedyMatching();

	runEdmonds(0);

      } else runEdmonds(1);

    }

    ///Runs Edmonds' algorithm.

    ///If heur=0 it runs Edmonds' algorithm. If heur=1 it runs

    ///Edmonds' algorithm with a heuristic of postponing shrinks,

    ///giving a faster algorithm for dense graphs.  

    void runEdmonds( int heur = 1 ) {

      //each vertex is put to C

      for(NodeIt v(g); v!=INVALID; ++v)

	position.set(v,C);      

      typename Graph::template NodeMap<Node> ear(g,INVALID); 

      //undefined for the base nodes of the blossoms (i.e. for the

      //representative elements of UFE blossom) and for the nodes in C 

      UFECrossRef blossom_base(g);

      UFE blossom(blossom_base);

      UFECrossRef tree_base(g);

      UFE tree(tree_base);

      //If these UFE's would be members of the class then also

      //blossom_base and tree_base should be a member.

      //We build only one tree and the other vertices uncovered by the

      //matching belong to C. (They can be considered as singleton

      //trees.) If this tree can be augmented or no more

      //grow/augmentation/shrink is possible then we return to this

      //"for" cycle.

      for(NodeIt v(g); v!=INVALID; ++v) {

	if ( position[v]==C && _mate[v]==INVALID ) {

	  blossom.insert(v);

	  tree.insert(v); 

	  position.set(v,D);

	  if ( heur == 1 ) lateShrink( v, ear, blossom, tree );

	  else normShrink( v, ear, blossom, tree );

	}

      }

    }

    ///Finds a greedy matching starting from the actual matching.

    ///Starting form the actual matching stored, it finds a maximal

    ///greedy matching.

    void greedyMatching() {

      for(NodeIt v(g); v!=INVALID; ++v)

	if ( _mate[v]==INVALID ) {

	  for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {

	    Node y=g.runningNode(e);

	    if ( _mate[y]==INVALID && y!=v ) {

	      _mate.set(v,y);

	      _mate.set(y,v);

	      break;

	    }

	  }

	} 

    }

    ///Returns the size of the actual matching stored.

    ///Returns the size of the actual matching stored. After \ref

    ///run() it returns the size of a maximum matching in the graph.

    int size() const {

      int s=0;

      for(NodeIt v(g); v!=INVALID; ++v) {

	if ( _mate[v]!=INVALID ) {

	  ++s;

	}

      }

      return s/2;

    }

    ///Resets the actual matching to the empty matching.

    ///Resets the actual matching to the empty matching.  

    ///

    void resetMatching() {

      for(NodeIt v(g); v!=INVALID; ++v)

	_mate.set(v,INVALID);      

    }

    ///Returns the mate of a node in the actual matching.

    ///Returns the mate of a \c node in the actual matching. 

    ///Returns INVALID if the \c node is not covered by the actual matching. 

    Node mate(Node& node) const {

      return _mate[node];

    } 

    ///Reads a matching from a \c Node valued \c Node map.

    ///Reads a matching from a \c Node valued \c Node map. This map

    ///must be \e symmetric, i.e. if \c map[u]==v then \c map[v]==u

    ///must hold, and \c uv will be an edge of the matching.

    template<typename NMapN>

    void readNMapNode(NMapN& map) {

      for(NodeIt v(g); v!=INVALID; ++v) {

	_mate.set(v,map[v]);   

      } 

    } 

    ///Writes the stored matching to a \c Node valued \c Node map.

    ///Writes the stored matching to a \c Node valued \c Node map. The

    ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c

    ///map[v]==u will hold, and now \c uv is an edge of the matching.

    template<typename NMapN>

    void writeNMapNode (NMapN& map) const {

      for(NodeIt v(g); v!=INVALID; ++v) {

	map.set(v,_mate[v]);   

      } 

    } 

    ///Reads a matching from an \c UEdge valued \c Node map.

    ///Reads a matching from an \c UEdge valued \c Node map. \c

    ///map[v] must be an \c UEdge incident to \c v. This map must

    ///have the property that if \c g.oppositeNode(u,map[u])==v then

    ///\c \c g.oppositeNode(v,map[v])==u holds, and now some edge

    ///joining \c u to \c v will be an edge of the matching.

    template<typename NMapE>

    void readNMapEdge(NMapE& map) {

      for(NodeIt v(g); v!=INVALID; ++v) {

	UEdge e=map[v];

	if ( e!=INVALID )

	  _mate.set(v,g.oppositeNode(v,e));

      } 

    } 

    ///Writes the matching stored to an \c UEdge valued \c Node map.

    ///Writes the stored matching to an \c UEdge valued \c Node

    ///map. \c map[v] will be an \c UEdge incident to \c v. This

    ///map will have the property that if \c g.oppositeNode(u,map[u])

    ///== v then \c map[u]==map[v] holds, and now this edge is an edge

    ///of the matching.

    template<typename NMapE>

    void writeNMapEdge (NMapE& map)  const {

      typename Graph::template NodeMap<bool> todo(g,true); 

      for(NodeIt v(g); v!=INVALID; ++v) {

	if ( todo[v] && _mate[v]!=INVALID ) {

	  Node u=_mate[v];

	  for(IncEdgeIt e(g,v); e!=INVALID; ++e) {

	    if ( g.runningNode(e) == u ) {

	      map.set(u,e);

	      map.set(v,e);

	      todo.set(u,false);

	      todo.set(v,false);

	      break;

	    }

	  }

	}

      } 

    }

    ///Reads a matching from a \c bool valued \c Edge map.

    ///Reads a matching from a \c bool valued \c Edge map. This map

    ///must have the property that there are no two incident edges \c

    ///e, \c f with \c map[e]==map[f]==true. The edges \c e with \c

    ///map[e]==true form the matching.

    template<typename EMapB>

    void readEMapBool(EMapB& map) {

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

	if ( map[e] ) {

	  Node u=g.source(e);	  

	  Node v=g.target(e);

	  _mate.set(u,v);

	  _mate.set(v,u);

	} 

      } 

    }

    ///Writes the matching stored to a \c bool valued \c Edge map.

    ///Writes the matching stored to a \c bool valued \c Edge

    ///map. This map will have the property that there are no two

    ///incident edges \c e, \c f with \c map[e]==map[f]==true. The

    ///edges \c e with \c map[e]==true form the matching.

    template<typename EMapB>

    void writeEMapBool (EMapB& map) const {

      for(UEdgeIt e(g); e!=INVALID; ++e) map.set(e,false);

      typename Graph::template NodeMap<bool> todo(g,true); 

      for(NodeIt v(g); v!=INVALID; ++v) {

	if ( todo[v] && _mate[v]!=INVALID ) {

	  Node u=_mate[v];

	  for(IncEdgeIt e(g,v); e!=INVALID; ++e) {

	    if ( g.runningNode(e) == u ) {

	      map.set(e,true);

	      todo.set(u,false);

	      todo.set(v,false);

	      break;

	    }

	  }

	}

      } 

    }

    ///Writes the canonical decomposition of the graph after running

    ///the algorithm.

    ///After calling any run methods of the class, it writes the

    ///Gallai-Edmonds canonical decomposition of the graph. \c map

    ///must be a node map of \ref pos_enum 's.

    template<typename NMapEnum>

    void writePos (NMapEnum& map) const {

      for(NodeIt v(g); v!=INVALID; ++v)  map.set(v,position[v]);

    }

  private: 

    void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,  

		    UFE& blossom, UFE& tree);

    void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,  

		    UFE& blossom, UFE& tree);

    void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear,  

		UFE& blossom, UFE& tree,std::queue<Node>& Q);

    void shrinkStep(Node& top, Node& middle, Node& bottom,

		    typename Graph::template NodeMap<Node>& ear,  

		    UFE& blossom, UFE& tree, std::queue<Node>& Q);

    bool growOrAugment(Node& y, Node& x, typename Graph::template 

		       NodeMap<Node>& ear, UFE& blossom, UFE& tree, 

		       std::queue<Node>& Q);

    void augment(Node x, typename Graph::template NodeMap<Node>& ear,  

		 UFE& blossom, UFE& tree);

  };

  // **********************************************************************

  //  IMPLEMENTATIONS

  // **********************************************************************

  template <typename Graph>

  void MaxMatching<Graph>::lateShrink(Node v, typename Graph::template 

				      NodeMap<Node>& ear, UFE& blossom, UFE& tree) {

    //We have one tree which we grow, and also shrink but only if it cannot be

    //postponed. If we augment then we return to the "for" cycle of

    //runEdmonds().

    std::queue<Node> Q;   //queue of the totally unscanned nodes

    Q.push(v);  

    std::queue<Node> R;   

    //queue of the nodes which must be scanned for a possible shrink

    while ( !Q.empty() ) {

      Node x=Q.front();

      Q.pop();

      for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {

	Node y=g.runningNode(e);

	//growOrAugment grows if y is covered by the matching and

	//augments if not. In this latter case it returns 1.

	if ( position[y]==C && growOrAugment(y, x, ear, blossom, tree, Q) ) return;

      }

      R.push(x);

    }

    while ( !R.empty() ) {

      Node x=R.front();

      R.pop();

      for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {

	Node y=g.runningNode(e);

	if ( position[y] == D && blossom.find(x) != blossom.find(y) )  

	  //Recall that we have only one tree.

	  shrink( x, y, ear, blossom, tree, Q);	

	while ( !Q.empty() ) {

	  Node x=Q.front();

	  Q.pop();

	  for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {

	    Node y=g.runningNode(e);

	    //growOrAugment grows if y is covered by the matching and

	    //augments if not. In this latter case it returns 1.

	    if ( position[y]==C && growOrAugment(y, x, ear, blossom, tree, Q) ) return;

	  }

	  R.push(x);

	}

      } //for e

    } // while ( !R.empty() )

  }

  template <typename Graph>

  void MaxMatching<Graph>::normShrink(Node v,

				      typename Graph::template

				      NodeMap<Node>& ear,  

				      UFE& blossom, UFE& tree) {

    //We have one tree, which we grow and shrink. If we augment then we

    //return to the "for" cycle of runEdmonds().

    std::queue<Node> Q;   //queue of the unscanned nodes

    Q.push(v);  

    while ( !Q.empty() ) {

      Node x=Q.front();

      Q.pop();

      for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {

	Node y=g.runningNode(e);

	switch ( position[y] ) {

	case D:          //x and y must be in the same tree

	  if ( blossom.find(x) != blossom.find(y) )

	    //x and y are in the same tree

	    shrink( x, y, ear, blossom, tree, Q);

	  break;

	case C:

	  //growOrAugment grows if y is covered by the matching and

	  //augments if not. In this latter case it returns 1.

	  if ( growOrAugment(y, x, ear, blossom, tree, Q) ) return;

	  break;

	default: break;

       	}

      }

    }

  }

  template <typename Graph>

    void MaxMatching<Graph>::shrink(Node x,Node y, typename 

				    Graph::template NodeMap<Node>& ear,  

				    UFE& blossom, UFE& tree, std::queue<Node>& Q) {

    //x and y are the two adjacent vertices in two blossoms.

    typename Graph::template NodeMap<bool> path(g,false);

    Node b=blossom.find(x);

    path.set(b,true);

    b=_mate[b];

    while ( b!=INVALID ) { 

      b=blossom.find(ear[b]);

      path.set(b,true);

      b=_mate[b];

    } //we go until the root through bases of blossoms and odd vertices

    Node top=y;

    Node middle=blossom.find(top);

    Node bottom=x;

    while ( !path[middle] )

      shrinkStep(top, middle, bottom, ear, blossom, tree, Q);

    //Until we arrive to a node on the path, we update blossom, tree

    //and the positions of the odd nodes.

    Node base=middle;

    top=x;

    middle=blossom.find(top);

    bottom=y;

    Node blossom_base=blossom.find(base);

    while ( middle!=blossom_base )

      shrinkStep(top, middle, bottom, ear, blossom, tree, Q);

    //Until we arrive to a node on the path, we update blossom, tree

    //and the positions of the odd nodes.

    blossom.makeRep(base);

  }

  template <typename Graph>

  void MaxMatching<Graph>::shrinkStep(Node& top, Node& middle, Node& bottom,

				      typename Graph::template

				      NodeMap<Node>& ear,  

				      UFE& blossom, UFE& tree,

				      std::queue<Node>& Q) {

    //We traverse a blossom and update everything.

    ear.set(top,bottom);

    Node t=top;

    while ( t!=middle ) {

      Node u=_mate[t];

      t=ear[u];

      ear.set(t,u);

    } 

    bottom=_mate[middle];

    position.set(bottom,D);

    Q.push(bottom);

    top=ear[bottom];		

    Node oldmiddle=middle;

    middle=blossom.find(top);

    tree.erase(bottom);

    tree.erase(oldmiddle);

    blossom.insert(bottom);

    blossom.join(bottom, oldmiddle);

    blossom.join(top, oldmiddle);

  }

  template <typename Graph>

  bool MaxMatching<Graph>::growOrAugment(Node& y, Node& x, typename Graph::template

					 NodeMap<Node>& ear, UFE& blossom, UFE& tree,

					 std::queue<Node>& Q) {

    //x is in a blossom in the tree, y is outside. If y is covered by

    //the matching we grow, otherwise we augment. In this case we

    //return 1.

    if ( _mate[y]!=INVALID ) {       //grow

      ear.set(y,x);

      Node w=_mate[y];

      blossom.insert(w);

      position.set(y,A);

      position.set(w,D);

      tree.insert(y);

      tree.insert(w);

      tree.join(y,blossom.find(x));  

      tree.join(w,y);  

      Q.push(w);

    } else {                      //augment 

      augment(x, ear, blossom, tree);

      _mate.set(x,y);

      _mate.set(y,x);

      return true;

    }

    return false;

  }

  template <typename Graph>

  void MaxMatching<Graph>::augment(Node x,

				   typename Graph::template NodeMap<Node>& ear,  

				   UFE& blossom, UFE& tree) { 

    Node v=_mate[x];

    while ( v!=INVALID ) {

      Node u=ear[v];

      _mate.set(v,u);

      Node tmp=v;

      v=_mate[u];

      _mate.set(u,tmp);

    }

    Node y=blossom.find(x);

    for (typename UFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {   

      if ( position[tit] == D ) {

	for (typename UFE::ItemIt bit(blossom, tit); bit != INVALID; ++bit) {  

	  position.set( bit ,C);

	}

	blossom.eraseClass(tit);

      } else position.set( tit ,C);

    }

    tree.eraseClass(y);

  }

} //END OF NAMESPACE LEMON

#endif //LEMON_MAX_MATCHING_H