source: lemon-0.x/lemon/pr_bipartite_matching.h

/* -*- 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);
    bpm.run();
    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);
    bpm.run();
    return bpm.checkedRunPerfect();
  }

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