1.1 --- a/lemon/Makefile.am	Thu Jan 25 14:36:21 2007 +0000
     1.2 +++ b/lemon/Makefile.am	Thu Jan 25 14:38:55 2007 +0000
     1.3 @@ -37,6 +37,7 @@
     1.4  	lemon/bfs.h \
     1.5  	lemon/bin_heap.h \
     1.6  	lemon/bipartite_matching.h \
     1.7 +	lemon/bp_matching.h \
     1.8  	lemon/bpugraph_adaptor.h \
     1.9  	lemon/bucket_heap.h \
    1.10  	lemon/color.h \

     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/lemon/bp_matching.h	Thu Jan 25 14:38:55 2007 +0000
     2.3 @@ -0,0 +1,381 @@
     2.4 +/* -*- C++ -*-
     2.5 + * lemon/preflow_matching.h - Part of LEMON, a generic C++ optimization library
     2.6 + *
     2.7 + * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     2.8 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
     2.9 + *
    2.10 + * Permission to use, modify and distribute this software is granted
    2.11 + * provided that this copyright notice appears in all copies. For
    2.12 + * precise terms see the accompanying LICENSE file.
    2.13 + *
    2.14 + * This software is provided "AS IS" with no warranty of any kind,
    2.15 + * express or implied, and with no claim as to its suitability for any
    2.16 + * purpose.
    2.17 + *
    2.18 + */
    2.19 +
    2.20 +#ifndef LEMON_BP_MATCHING
    2.21 +#define LEMON_BP_MATCHING
    2.22 +
    2.23 +#include <lemon/graph_utils.h>
    2.24 +#include <lemon/iterable_maps.h>
    2.25 +#include <iostream>
    2.26 +#include <queue>
    2.27 +#include <lemon/counter.h>
    2.28 +#include <lemon/elevator.h>
    2.29 +
    2.30 +///\ingroup matching
    2.31 +///\file
    2.32 +///\brief Push-prelabel maximum matching algorithms in bipartite graphs.
    2.33 +///
    2.34 +///\todo This file slightly conflicts with \ref lemon/bipartite_matching.h
    2.35 +///\todo (Re)move the XYZ_TYPEDEFS macros
    2.36 +namespace lemon {
    2.37 +
    2.38 +#define BIPARTITE_TYPEDEFS(Graph)		\
    2.39 +  GRAPH_TYPEDEFS(Graph)				\
    2.40 +    typedef Graph::ANodeIt ANodeIt;	\
    2.41 +    typedef Graph::BNodeIt BNodeIt;
    2.42 +
    2.43 +#define UNDIRBIPARTITE_TYPEDEFS(Graph)		\
    2.44 +  UNDIRGRAPH_TYPEDEFS(Graph)			\
    2.45 +    typedef Graph::ANodeIt ANodeIt;	\
    2.46 +    typedef Graph::BNodeIt BNodeIt;
    2.47 +
    2.48 +  template<class Graph,
    2.49 +	   class MT=typename Graph::template ANodeMap<typename Graph::UEdge> >
    2.50 +  class BpMatching {
    2.51 +    typedef typename Graph::Node Node;
    2.52 +    typedef typename Graph::ANodeIt ANodeIt;
    2.53 +    typedef typename Graph::BNodeIt BNodeIt;
    2.54 +    typedef typename Graph::UEdge UEdge;
    2.55 +    typedef typename Graph::IncEdgeIt IncEdgeIt;
    2.56 +    
    2.57 +    const Graph &_g;
    2.58 +    int _node_num;
    2.59 +    MT &_matching;
    2.60 +    Elevator<Graph,typename Graph::BNode> _levels;
    2.61 +    typename Graph::template BNodeMap<int> _cov;
    2.62 +
    2.63 +  public:
    2.64 +    BpMatching(const Graph &g, MT &matching) :
    2.65 +      _g(g),
    2.66 +      _node_num(countBNodes(g)),
    2.67 +      _matching(matching),
    2.68 +      _levels(g,_node_num),
    2.69 +      _cov(g,0)
    2.70 +    {
    2.71 +    }
    2.72 +    
    2.73 +  private:
    2.74 +    void init() 
    2.75 +    {
    2.76 +//     for(BNodeIt n(g);n!=INVALID;++n) cov[n]=0;
    2.77 +      for(ANodeIt n(_g);n!=INVALID;++n)
    2.78 +	if((_matching[n]=IncEdgeIt(_g,n))!=INVALID)
    2.79 +	  ++_cov[_g.oppositeNode(n,_matching[n])];
    2.80 +
    2.81 +      std::queue<Node> q;
    2.82 +      _levels.initStart();
    2.83 +      for(BNodeIt n(_g);n!=INVALID;++n)
    2.84 +	if(_cov[n]>1) {
    2.85 +	  _levels.initAddItem(n);
    2.86 +	  q.push(n);
    2.87 +	}
    2.88 +      int hlev=0;
    2.89 +      while(!q.empty()) {
    2.90 +	Node n=q.front();
    2.91 +	q.pop();
    2.92 +	int nlev=_levels[n]+1;
    2.93 +	for(IncEdgeIt e(_g,n);e!=INVALID;++e) {
    2.94 +	  Node m=_g.runningNode(e);
    2.95 +	  if(e==_matching[m]) {
    2.96 +	    for(IncEdgeIt f(_g,m);f!=INVALID;++f) {
    2.97 +	      Node r=_g.runningNode(f);
    2.98 +	      if(_levels[r]>nlev) {
    2.99 +		for(;nlev>hlev;hlev++)
   2.100 +		  _levels.initNewLevel();
   2.101 +		_levels.initAddItem(r);
   2.102 +		q.push(r);
   2.103 +	      }
   2.104 +	    }
   2.105 +	  }
   2.106 +	}
   2.107 +      }
   2.108 +      _levels.initFinish();
   2.109 +      for(BNodeIt n(_g);n!=INVALID;++n)
   2.110 +	if(_cov[n]<1&&_levels[n]<_node_num)
   2.111 +	  _levels.activate(n);
   2.112 +    }
   2.113 +  public:
   2.114 +    int run() 
   2.115 +    {
   2.116 +      init();
   2.117 +
   2.118 +      Node act;
   2.119 +      Node bact=INVALID;
   2.120 +      Node last_activated=INVALID;
   2.121 +//       while((act=last_activated!=INVALID?
   2.122 +// 	     last_activated:_levels.highestActive())
   2.123 +// 	    !=INVALID)
   2.124 +      while((act=_levels.highestActive())!=INVALID) {
   2.125 +	last_activated=INVALID;
   2.126 +	int actlevel=_levels[act];
   2.127 +	
   2.128 +	UEdge bedge=INVALID;
   2.129 +	int nlevel=_node_num;
   2.130 +	{
   2.131 +	  int nnlevel;
   2.132 +	  for(IncEdgeIt tbedge(_g,act);
   2.133 +	      tbedge!=INVALID && nlevel>=actlevel;
   2.134 +	      ++tbedge)
   2.135 +	    if((nnlevel=_levels[_g.bNode(_matching[_g.runningNode(tbedge)])])<
   2.136 +	       nlevel)
   2.137 +	      {
   2.138 +		nlevel=nnlevel;
   2.139 +		bedge=tbedge;
   2.140 +	      }
   2.141 +	}
   2.142 +	if(nlevel<_node_num) {
   2.143 +	  if(nlevel>=actlevel)
   2.144 +	    _levels.liftHighestActiveTo(nlevel+1);
   2.145 +// 	    _levels.liftTo(act,nlevel+1);
   2.146 +	  bact=_g.bNode(_matching[_g.aNode(bedge)]);
   2.147 +	  if(--_cov[bact]<1) {
   2.148 +	    _levels.activate(bact);
   2.149 +	    last_activated=bact;
   2.150 +	  }
   2.151 +	  _matching[_g.aNode(bedge)]=bedge;
   2.152 +	  _cov[act]=1;
   2.153 +	  _levels.deactivate(act);
   2.154 +	}
   2.155 +	else {
   2.156 +	  if(_node_num>actlevel) 
   2.157 +	    _levels.liftHighestActiveTo(_node_num);
   2.158 +//  	    _levels.liftTo(act,_node_num);
   2.159 +	  _levels.deactivate(act); 
   2.160 +	}
   2.161 +
   2.162 +	if(_levels.onLevel(actlevel)==0)
   2.163 +	  _levels.liftToTop(actlevel);
   2.164 +      }
   2.165 +      
   2.166 +      int ret=_node_num;
   2.167 +      for(ANodeIt n(_g);n!=INVALID;++n)
   2.168 +	if(_matching[n]==INVALID) ret--;
   2.169 +	else if (_cov[_g.bNode(_matching[n])]>1) {
   2.170 +	  _cov[_g.bNode(_matching[n])]--;
   2.171 +	  ret--;
   2.172 +	  _matching[n]=INVALID;
   2.173 +	}
   2.174 +      return ret;
   2.175 +    }
   2.176 +    
   2.177 +    ///\returns -1 if there is a perfect matching, or an empty level
   2.178 +    ///if it doesn't exists
   2.179 +    int runPerfect() 
   2.180 +    {
   2.181 +      init();
   2.182 +
   2.183 +      Node act;
   2.184 +      Node bact=INVALID;
   2.185 +      Node last_activated=INVALID;
   2.186 +      while((act=_levels.highestActive())!=INVALID) {
   2.187 +	last_activated=INVALID;
   2.188 +	int actlevel=_levels[act];
   2.189 +	
   2.190 +	UEdge bedge=INVALID;
   2.191 +	int nlevel=_node_num;
   2.192 +	{
   2.193 +	  int nnlevel;
   2.194 +	  for(IncEdgeIt tbedge(_g,act);
   2.195 +	      tbedge!=INVALID && nlevel>=actlevel;
   2.196 +	      ++tbedge)
   2.197 +	    if((nnlevel=_levels[_g.bNode(_matching[_g.runningNode(tbedge)])])<
   2.198 +	       nlevel)
   2.199 +	      {
   2.200 +		nlevel=nnlevel;
   2.201 +		bedge=tbedge;
   2.202 +	      }
   2.203 +	}
   2.204 +	if(nlevel<_node_num) {
   2.205 +	  if(nlevel>=actlevel)
   2.206 +	    _levels.liftHighestActiveTo(nlevel+1);
   2.207 +	  bact=_g.bNode(_matching[_g.aNode(bedge)]);
   2.208 +	  if(--_cov[bact]<1) {
   2.209 +	    _levels.activate(bact);
   2.210 +	    last_activated=bact;
   2.211 +	  }
   2.212 +	  _matching[_g.aNode(bedge)]=bedge;
   2.213 +	  _cov[act]=1;
   2.214 +	  _levels.deactivate(act);
   2.215 +	}
   2.216 +	else {
   2.217 +	  if(_node_num>actlevel) 
   2.218 +	    _levels.liftHighestActiveTo(_node_num);
   2.219 +	  _levels.deactivate(act); 
   2.220 +	}
   2.221 +
   2.222 +	if(_levels.onLevel(actlevel)==0)
   2.223 +	  return actlevel;
   2.224 +      }
   2.225 +      return -1;
   2.226 +    }
   2.227 + 
   2.228 +    template<class GT>
   2.229 +    void aBarrier(GT &bar,int empty_level=-1) 
   2.230 +    {
   2.231 +      if(empty_level==-1)
   2.232 +	for(empty_level=0;_levels.onLevel(empty_level);empty_level++) ;
   2.233 +      for(ANodeIt n(_g);n!=INVALID;++n)
   2.234 +	bar[n] = _matching[n]==INVALID ||
   2.235 +	  _levels[_g.bNode(_matching[n])]<empty_level;  
   2.236 +    }  
   2.237 +    template<class GT>
   2.238 +    void bBarrier(GT &bar, int empty_level=-1) 
   2.239 +    {
   2.240 +      if(empty_level==-1)
   2.241 +	for(empty_level=0;_levels.onLevel(empty_level);empty_level++) ;
   2.242 +      for(BNodeIt n(_g);n!=INVALID;++n) bar[n]=(_levels[n]>empty_level);  
   2.243 +    }  
   2.244 +  
   2.245 +  };
   2.246 +  
   2.247 +  
   2.248 +  ///Maximum cardinality of the matchings in a bipartite graph
   2.249 +
   2.250 +  ///\ingroup matching
   2.251 +  ///This function finds the maximum cardinality of the matchings
   2.252 +  ///in a bipartite graph \c g.
   2.253 +  ///\param g An undirected bipartite graph.
   2.254 +  ///\return The cardinality of the maximum matching.
   2.255 +  ///
   2.256 +  ///\note The the implementation is based
   2.257 +  ///on the push-relabel principle.
   2.258 +  template<class Graph>
   2.259 +  int maxBpMatching(const Graph &g)
   2.260 +  {
   2.261 +    typename Graph::template ANodeMap<typename Graph::UEdge> matching(g);
   2.262 +    return maxBpMatching(g,matching);
   2.263 +  }
   2.264 +
   2.265 +  ///Maximum cardinality matching in a bipartite graph
   2.266 +
   2.267 +  ///\ingroup matching
   2.268 +  ///This function finds a maximum cardinality matching
   2.269 +  ///in a bipartite graph \c g.
   2.270 +  ///\param g An undirected bipartite graph.
   2.271 +  ///\retval matching A readwrite ANodeMap of value type \c Edge.
   2.272 +  /// The found edges will be returned in this map,
   2.273 +  /// i.e. for an \c ANode \c n,
   2.274 +  /// the edge <tt>matching[n]</tt> is the one that covers the node \c n, or
   2.275 +  /// \ref INVALID if it is uncovered.
   2.276 +  ///\return The cardinality of the maximum matching.
   2.277 +  ///
   2.278 +  ///\note The the implementation is based
   2.279 +  ///on the push-relabel principle.
   2.280 +  template<class Graph,class MT>
   2.281 +  int maxBpMatching(const Graph &g,MT &matching) 
   2.282 +  {
   2.283 +    return BpMatching<Graph,MT>(g,matching).run();
   2.284 +  }
   2.285 +
   2.286 +  ///Maximum cardinality matching in a bipartite graph
   2.287 +
   2.288 +  ///\ingroup matching
   2.289 +  ///This function finds a maximum cardinality matching
   2.290 +  ///in a bipartite graph \c g.
   2.291 +  ///\param g An undirected bipartite graph.
   2.292 +  ///\retval matching A readwrite ANodeMap of value type \c Edge.
   2.293 +  /// The found edges will be returned in this map,
   2.294 +  /// i.e. for an \c ANode \c n,
   2.295 +  /// the edge <tt>matching[n]</tt> is the one that covers the node \c n, or
   2.296 +  /// \ref INVALID if it is uncovered.
   2.297 +  ///\retval barrier A \c bool WriteMap on the BNodes. The map will be set
   2.298 +  /// exactly once for each BNode. The nodes with \c true value represent
   2.299 +  /// a barrier \e B, i.e. the cardinality of \e B minus the number of its
   2.300 +  /// neighbor is equal to the number of the <tt>BNode</tt>s minus the
   2.301 +  /// cardinality of the maximum matching.
   2.302 +  ///\return The cardinality of the maximum matching.
   2.303 +  ///
   2.304 +  ///\note The the implementation is based
   2.305 +  ///on the push-relabel principle.
   2.306 +  template<class Graph,class MT, class GT>
   2.307 +  int maxBpMatching(const Graph &g,MT &matching,GT &barrier) 
   2.308 +  {
   2.309 +    BpMatching<Graph,MT> bpm(g,matching);
   2.310 +    int ret=bpm.run();
   2.311 +    bpm.barrier(barrier);
   2.312 +    return ret;
   2.313 +  }  
   2.314 +
   2.315 +  ///Perfect matching in a bipartite graph
   2.316 +
   2.317 +  ///\ingroup matching
   2.318 +  ///This function checks whether the bipartite graph \c g
   2.319 +  ///has a perfect matching.
   2.320 +  ///\param g An undirected bipartite graph.
   2.321 +  ///\return \c true iff \c g has a perfect matching.
   2.322 +  ///
   2.323 +  ///\note The the implementation is based
   2.324 +  ///on the push-relabel principle.
   2.325 +  template<class Graph>
   2.326 +  bool perfectBpMatching(const Graph &g)
   2.327 +  {
   2.328 +    typename Graph::template ANodeMap<typename Graph::UEdge> matching(g);
   2.329 +    return perfectBpMatching(g,matching);
   2.330 +  }
   2.331 +
   2.332 +  ///Perfect matching in a bipartite graph
   2.333 +
   2.334 +  ///\ingroup matching
   2.335 +  ///This function finds a perfect matching in a bipartite graph \c g.
   2.336 +  ///\param g An undirected bipartite graph.
   2.337 +  ///\retval matching A readwrite ANodeMap of value type \c Edge.
   2.338 +  /// The found edges will be returned in this map,
   2.339 +  /// i.e. for an \c ANode \c n,
   2.340 +  /// the edge <tt>matching[n]</tt> is the one that covers the node \c n.
   2.341 +  /// The values are unspecified if the graph
   2.342 +  /// has no perfect matching.
   2.343 +  ///\return \c true iff \c g has a perfect matching.
   2.344 +  ///
   2.345 +  ///\note The the implementation is based
   2.346 +  ///on the push-relabel principle.
   2.347 +  template<class Graph,class MT>
   2.348 +  bool perfectBpMatching(const Graph &g,MT &matching) 
   2.349 +  {
   2.350 +    return BpMatching<Graph,MT>(g,matching).runPerfect()<0;
   2.351 +  }
   2.352 +
   2.353 +  ///Perfect matching in a bipartite graph
   2.354 +
   2.355 +  ///\ingroup matching
   2.356 +  ///This function finds a perfect matching in a bipartite graph \c g.
   2.357 +  ///\param g An undirected bipartite graph.
   2.358 +  ///\retval matching A readwrite ANodeMap of value type \c Edge.
   2.359 +  /// The found edges will be returned in this map,
   2.360 +  /// i.e. for an \c ANode \c n,
   2.361 +  /// the edge <tt>matching[n]</tt> is the one that covers the node \c n.
   2.362 +  /// The values are unspecified if the graph
   2.363 +  /// has no perfect matching.
   2.364 +  ///\retval barrier A \c bool WriteMap on the BNodes. The map will only
   2.365 +  /// be set if \c g has no perfect matching. In this case it is set 
   2.366 +  /// exactly once for each BNode. The nodes with \c true value represent
   2.367 +  /// a barrier, i.e. a subset \e B a of BNodes with the property that
   2.368 +  /// the cardinality of \e B is greater than the numner of its neighbors.
   2.369 +  ///\return \c true iff \c g has a perfect matching.
   2.370 +  ///
   2.371 +  ///\note The the implementation is based
   2.372 +  ///on the push-relabel principle.
   2.373 +  template<class Graph,class MT, class GT>
   2.374 +  int perfectBpMatching(const Graph &g,MT &matching,GT &barrier) 
   2.375 +  {
   2.376 +    BpMatching<Graph,MT> bpm(g,matching);
   2.377 +    int ret=bpm.run();
   2.378 +    if(ret>=0)
   2.379 +      bpm.barrier(barrier,ret);
   2.380 +    return ret<0;
   2.381 +  }  
   2.382 +}
   2.383 +
   2.384 +#endif