diff -r 1dd4d6ff9bac -r 7a096a6bf53a lemon/pr_bipartite_matching.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/pr_bipartite_matching.h	Sat Aug 11 16:34:41 2007 +0000
@@ -0,0 +1,572 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * 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:
+
+    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.liftHighestActiveTo(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 {
+	  if(_node_num>actlevel) 
+	    _levels.liftHighestActiveTo(_node_num);
+	  _levels.deactivate(act); 
+	}
+
+	if(_levels.onLevel(actlevel)==0)
+	  _levels.liftToTop(actlevel);
+      }
+      
+      _matching_size = _node_num;
+      for(ANodeIt n(_g);n!=INVALID;++n)
+	if(_matching[n]==INVALID) _matching_size--;
+	else if (_cov[_g.bNode(_matching[n])]>1) {
+	  _cov[_g.bNode(_matching[n])]--;
+	  _matching_size--;
+	  _matching[n]=INVALID;
+	}
+    }
+
+    ///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.liftHighestActiveTo(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 {
+	  if(_node_num>actlevel) 
+	    _levels.liftHighestActiveTo(_node_num);
+	  _levels.deactivate(act); 
+	}
+
+	if(_levels.onLevel(actlevel)==0)
+	  _empty_level=actlevel;
+	  return false;
+      }
+      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.
+    ///@{
+
+    /// \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& mm) const {
+      for (ANodeIt n(_g);n!=INVALID;++n) {
+        if (_matching[n]!=INVALID) mm.set(_matching[n],true);
+      }
+      return _matching_size;
+    }
+
+    ///\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<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;
+    }
+
+
+    ///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 UEdgeMap of value type \c bool.
+  /// The found edges will be returned in this map.
+  ///\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.matching(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 UEdgeMap of value type \c bool.
+  /// The found edges will be returned in this map.
+  ///\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.matching(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 UEdgeMap of value type \c bool.
+  /// The found edges will be returned in this map.
+  /// 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.runPerfect();
+    if (ret) bpm.matching(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 readwrite UEdgeMap of value type \c bool.
+  /// The found edges will be returned in this map.
+  /// 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>
+  int prPerfectBipartiteMatching(const Graph &g,MT &matching,GT &barrier) 
+  {
+    PrBipartiteMatching<Graph> bpm(g);
+    bool ret=bpm.runPerfect();
+    if(ret)
+      bpm.matching(matching);
+    else
+      bpm.bBarrier(barrier);
+    return ret;
+  }  
+}
+
+#endif