src/work/edmonds_karp.h
author marci
Wed, 17 Mar 2004 18:18:26 +0000
changeset 198 5cec393baade
parent 197 fff43d9c7110
child 206 47f62d789fe7
permissions -rw-r--r--
max cardinality bipartite matching demo, something to play with it
     1 // -*- c++ -*-
     2 #ifndef EDMONDS_KARP_H
     3 #define EDMONDS_KARP_H
     4 
     5 #include <algorithm>
     6 #include <list>
     7 #include <iterator>
     8 
     9 #include <bfs_iterator.h>
    10 #include <invalid.h>
    11 
    12 namespace hugo {
    13 
    14   template<typename Graph, typename Number, typename FlowMap, typename CapacityMap>
    15   class ResGraph {
    16   public:
    17     typedef typename Graph::Node Node;
    18     typedef typename Graph::NodeIt NodeIt;
    19   private:
    20     typedef typename Graph::SymEdgeIt OldSymEdgeIt;
    21     const Graph& G;
    22     FlowMap& flow;
    23     const CapacityMap& capacity;
    24   public:
    25     ResGraph(const Graph& _G, FlowMap& _flow, 
    26 	     const CapacityMap& _capacity) : 
    27       G(_G), flow(_flow), capacity(_capacity) { }
    28 
    29     class Edge; 
    30     class OutEdgeIt; 
    31     friend class Edge; 
    32     friend class OutEdgeIt; 
    33 
    34     class Edge {
    35       friend class ResGraph<Graph, Number, FlowMap, CapacityMap>;
    36     protected:
    37       const ResGraph<Graph, Number, FlowMap, CapacityMap>* resG;
    38       OldSymEdgeIt sym;
    39     public:
    40       Edge() { } 
    41       //Edge(const Edge& e) : resG(e.resG), sym(e.sym) { }
    42       Number free() const { 
    43 	if (resG->G.aNode(sym)==resG->G.tail(sym)) { 
    44 	  return (resG->capacity.get(sym)-resG->flow.get(sym)); 
    45 	} else { 
    46 	  return (resG->flow.get(sym)); 
    47 	}
    48       }
    49       bool valid() const { return sym.valid(); }
    50       void augment(Number a) const {
    51 	if (resG->G.aNode(sym)==resG->G.tail(sym)) { 
    52 	  resG->flow.set(sym, resG->flow.get(sym)+a);
    53 	  //resG->flow[sym]+=a;
    54 	} else { 
    55 	  resG->flow.set(sym, resG->flow.get(sym)-a);
    56 	  //resG->flow[sym]-=a;
    57 	}
    58       }
    59     };
    60 
    61     class OutEdgeIt : public Edge {
    62       friend class ResGraph<Graph, Number, FlowMap, CapacityMap>;
    63     public:
    64       OutEdgeIt() { }
    65       //OutEdgeIt(const OutEdgeIt& e) { resG=e.resG; sym=e.sym; }
    66     private:
    67       OutEdgeIt(const ResGraph<Graph, Number, FlowMap, CapacityMap>& _resG, Node v) { 
    68       	resG=&_resG;
    69 	sym=resG->G.template first<OldSymEdgeIt>(v);
    70 	while( sym.valid() && !(free()>0) ) { ++sym; }
    71       }
    72     public:
    73       OutEdgeIt& operator++() { 
    74 	++sym; 
    75 	while( sym.valid() && !(free()>0) ) { ++sym; }
    76 	return *this; 
    77       }
    78     };
    79 
    80     void /*getF*/first(OutEdgeIt& e, Node v) const { 
    81       e=OutEdgeIt(*this, v); 
    82     }
    83     void /*getF*/first(NodeIt& v) const { G./*getF*/first(v); }
    84     
    85     template< typename It >
    86     It first() const { 
    87       It e;      
    88       /*getF*/first(e);
    89       return e; 
    90     }
    91 
    92     template< typename It >
    93     It first(Node v) const { 
    94       It e;
    95       /*getF*/first(e, v);
    96       return e; 
    97     }
    98 
    99     Node tail(Edge e) const { return G.aNode(e.sym); }
   100     Node head(Edge e) const { return G.bNode(e.sym); }
   101 
   102     Node aNode(OutEdgeIt e) const { return G.aNode(e.sym); }
   103     Node bNode(OutEdgeIt e) const { return G.bNode(e.sym); }
   104 
   105     int id(Node v) const { return G.id(v); }
   106 
   107     template <typename S>
   108     class NodeMap {
   109       typename Graph::NodeMap<S> node_map; 
   110     public:
   111       NodeMap(const ResGraph<Graph, Number, FlowMap, CapacityMap>& _G) : node_map(_G.G) { }
   112       NodeMap(const ResGraph<Graph, Number, FlowMap, CapacityMap>& _G, S a) : node_map(_G.G, a) { }
   113       void set(Node nit, S a) { node_map.set(nit, a); }
   114       S get(Node nit) const { return node_map.get(nit); }
   115       S& operator[](Node nit) { return node_map[nit]; } 
   116       const S& operator[](Node nit) const { return node_map[nit]; } 
   117     };
   118 
   119   };
   120 
   121 
   122   template<typename Graph, typename Number, typename FlowMap, typename CapacityMap>
   123   class ResGraph2 {
   124   public:
   125     typedef typename Graph::Node Node;
   126     typedef typename Graph::NodeIt NodeIt;
   127   private:
   128     //typedef typename Graph::SymEdgeIt OldSymEdgeIt;
   129     typedef typename Graph::OutEdgeIt OldOutEdgeIt;
   130     typedef typename Graph::InEdgeIt OldInEdgeIt;
   131     
   132     const Graph& G;
   133     FlowMap& flow;
   134     const CapacityMap& capacity;
   135   public:
   136     ResGraph2(const Graph& _G, FlowMap& _flow, 
   137 	     const CapacityMap& _capacity) : 
   138       G(_G), flow(_flow), capacity(_capacity) { }
   139 
   140     class Edge; 
   141     class OutEdgeIt; 
   142     friend class Edge; 
   143     friend class OutEdgeIt; 
   144 
   145     class Edge {
   146       friend class ResGraph2<Graph, Number, FlowMap, CapacityMap>;
   147     protected:
   148       const ResGraph2<Graph, Number, FlowMap, CapacityMap>* resG;
   149       //OldSymEdgeIt sym;
   150       OldOutEdgeIt out;
   151       OldInEdgeIt in;
   152       bool out_or_in; //true, iff out
   153     public:
   154       Edge() : out_or_in(true) { } 
   155       Number free() const { 
   156 	if (out_or_in) { 
   157 	  return (resG->capacity.get(out)-resG->flow.get(out)); 
   158 	} else { 
   159 	  return (resG->flow.get(in)); 
   160 	}
   161       }
   162       bool valid() const { 
   163 	return out_or_in && out.valid() || in.valid(); }
   164       void augment(Number a) const {
   165 	if (out_or_in) { 
   166 	  resG->flow.set(out, resG->flow.get(out)+a);
   167 	} else { 
   168 	  resG->flow.set(in, resG->flow.get(in)-a);
   169 	}
   170       }
   171     };
   172 
   173     class OutEdgeIt : public Edge {
   174       friend class ResGraph2<Graph, Number, FlowMap, CapacityMap>;
   175     public:
   176       OutEdgeIt() { }
   177     private:
   178       OutEdgeIt(const ResGraph2<Graph, Number, FlowMap, CapacityMap>& _resG, Node v) { 
   179       	resG=&_resG;
   180 	out=resG->G.template first<OldOutEdgeIt>(v);
   181 	while( out.valid() && !(free()>0) ) { ++out; }
   182 	if (!out.valid()) {
   183 	  out_or_in=0;
   184 	  in=resG->G.template first<OldInEdgeIt>(v);
   185 	  while( in.valid() && !(free()>0) ) { ++in; }
   186 	}
   187       }
   188     public:
   189       OutEdgeIt& operator++() { 
   190 	if (out_or_in) {
   191 	  Node v=resG->G.aNode(out);
   192 	  ++out;
   193 	  while( out.valid() && !(free()>0) ) { ++out; }
   194 	  if (!out.valid()) {
   195 	    out_or_in=0;
   196 	    in=resG->G.template first<OldInEdgeIt>(v);
   197 	    while( in.valid() && !(free()>0) ) { ++in; }
   198 	  }
   199 	} else {
   200 	  ++in;
   201 	  while( in.valid() && !(free()>0) ) { ++in; } 
   202 	}
   203 	return *this; 
   204       }
   205     };
   206 
   207     void /*getF*/first(OutEdgeIt& e, Node v) const { 
   208       e=OutEdgeIt(*this, v); 
   209     }
   210     void /*getF*/first(NodeIt& v) const { G./*getF*/first(v); }
   211     
   212     template< typename It >
   213     It first() const { 
   214       It e;
   215       /*getF*/first(e);
   216       return e; 
   217     }
   218 
   219     template< typename It >
   220     It first(Node v) const { 
   221       It e;
   222       /*getF*/first(e, v);
   223       return e; 
   224     }
   225 
   226     Node tail(Edge e) const { 
   227       return ((e.out_or_in) ? G.aNode(e.out) : G.aNode(e.in)); }
   228     Node head(Edge e) const { 
   229       return ((e.out_or_in) ? G.bNode(e.out) : G.bNode(e.in)); }
   230 
   231     Node aNode(OutEdgeIt e) const { 
   232       return ((e.out_or_in) ? G.aNode(e.out) : G.aNode(e.in)); }
   233     Node bNode(OutEdgeIt e) const { 
   234       return ((e.out_or_in) ? G.bNode(e.out) : G.bNode(e.in)); }
   235 
   236     int id(Node v) const { return G.id(v); }
   237 
   238     template <typename S>
   239     class NodeMap {
   240       typename Graph::NodeMap<S> node_map; 
   241     public:
   242       NodeMap(const ResGraph2<Graph, Number, FlowMap, CapacityMap>& _G) : node_map(_G.G) { }
   243       NodeMap(const ResGraph2<Graph, Number, FlowMap, CapacityMap>& _G, S a) : node_map(_G.G, a) { }
   244       void set(Node nit, S a) { node_map.set(nit, a); }
   245       S get(Node nit) const { return node_map.get(nit); }
   246     };
   247   };
   248 
   249 
   250   template <typename Graph, typename Number, typename FlowMap, typename CapacityMap>
   251   class MaxFlow {
   252   public:
   253     typedef typename Graph::Node Node;
   254     typedef typename Graph::Edge Edge;
   255     typedef typename Graph::EdgeIt EdgeIt;
   256     typedef typename Graph::OutEdgeIt OutEdgeIt;
   257     typedef typename Graph::InEdgeIt InEdgeIt;
   258 
   259   private:
   260     const Graph* G;
   261     Node s;
   262     Node t;
   263     FlowMap* flow;
   264     const CapacityMap* capacity;
   265     typedef ResGraphWrapper<Graph, Number, FlowMap, CapacityMap > AugGraph;
   266     typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
   267     typedef typename AugGraph::Edge AugEdge;
   268 
   269   public:
   270     MaxFlow(const Graph& _G, Node _s, Node _t, FlowMap& _flow, const CapacityMap& _capacity) : 
   271       G(&_G), s(_s), t(_t), flow(&_flow), capacity(&_capacity) { }
   272     bool augmentOnShortestPath() {
   273       AugGraph res_graph(*G, *flow, *capacity);
   274       bool _augment=false;
   275       
   276       typedef typename AugGraph::NodeMap<bool> ReachedMap;
   277       BfsIterator5< AugGraph, /*AugOutEdgeIt,*/ ReachedMap > res_bfs(res_graph);
   278       res_bfs.pushAndSetReached(s);
   279 	
   280       typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
   281       pred.set(s, AugEdge(INVALID));
   282       
   283       typename AugGraph::NodeMap<Number> free(res_graph);
   284 	
   285       //searching for augmenting path
   286       while ( !res_bfs.finished() ) { 
   287 	AugOutEdgeIt e=res_bfs;
   288 	if (res_graph.valid(e) && res_bfs.isBNodeNewlyReached()) {
   289 	  Node v=res_graph.tail(e);
   290 	  Node w=res_graph.head(e);
   291 	  pred.set(w, e);
   292 	  if (res_graph.valid(pred.get(v))) {
   293 	    free.set(w, std::min(free.get(v), res_graph.free(e)));
   294 	  } else {
   295 	    free.set(w, res_graph.free(e)); 
   296 	  }
   297 	  if (res_graph.head(e)==t) { _augment=true; break; }
   298 	}
   299 	
   300 	++res_bfs;
   301       } //end of searching augmenting path
   302 
   303       if (_augment) {
   304 	Node n=t;
   305 	Number augment_value=free.get(t);
   306 	while (res_graph.valid(pred.get(n))) { 
   307 	  AugEdge e=pred.get(n);
   308 	  res_graph.augment(e, augment_value); 
   309 	  n=res_graph.tail(e);
   310 	}
   311       }
   312 
   313       return _augment;
   314     }
   315 
   316     template<typename MutableGraph> bool augmentOnBlockingFlow() {      
   317       bool _augment=false;
   318 
   319       AugGraph res_graph(*G, *flow, *capacity);
   320 
   321       typedef typename AugGraph::NodeMap<bool> ReachedMap;
   322       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > bfs(res_graph);
   323 
   324       bfs.pushAndSetReached(s);
   325       typename AugGraph::NodeMap<int> dist(res_graph); //filled up with 0's
   326       while ( !bfs.finished() ) { 
   327 	AugOutEdgeIt e=bfs;
   328 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   329 	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
   330 	}
   331 	
   332 	++bfs;
   333       } //computing distances from s in the residual graph
   334 
   335       MutableGraph F;
   336       typename AugGraph::NodeMap<typename MutableGraph::Node> 
   337 	res_graph_to_F(res_graph);
   338       for(typename AugGraph::NodeIt n=res_graph.template first<typename AugGraph::NodeIt>(); res_graph.valid(n); res_graph.next(n)) {
   339 	res_graph_to_F.set(n, F.addNode());
   340       }
   341       
   342       typename MutableGraph::Node sF=res_graph_to_F.get(s);
   343       typename MutableGraph::Node tF=res_graph_to_F.get(t);
   344 
   345       typename MutableGraph::EdgeMap<AugEdge> original_edge(F);
   346       typename MutableGraph::EdgeMap<Number> residual_capacity(F);
   347 
   348       //Making F to the graph containing the edges of the residual graph 
   349       //which are in some shortest paths
   350       for(typename AugGraph::EdgeIt e=res_graph.template first<typename AugGraph::EdgeIt>(); res_graph.valid(e); res_graph.next(e)) {
   351 	if (dist.get(res_graph.head(e))==dist.get(res_graph.tail(e))+1) {
   352 	  typename MutableGraph::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
   353 	  original_edge.update();
   354 	  original_edge.set(f, e);
   355 	  residual_capacity.update();
   356 	  residual_capacity.set(f, res_graph.free(e));
   357 	} 
   358       }
   359 
   360       bool __augment=true;
   361 
   362       while (__augment) {
   363 	__augment=false;
   364 	//computing blocking flow with dfs
   365 	typedef typename MutableGraph::NodeMap<bool> BlockingReachedMap;
   366 	DfsIterator4< MutableGraph, typename MutableGraph::OutEdgeIt, BlockingReachedMap > dfs(F);
   367 	typename MutableGraph::NodeMap<typename MutableGraph::Edge> pred(F);
   368 	pred.set(sF, typename MutableGraph::Edge(INVALID));
   369 	//invalid iterators for sources
   370 
   371 	typename MutableGraph::NodeMap<Number> free(F);
   372 
   373 	dfs.pushAndSetReached(sF);      
   374 	while (!dfs.finished()) {
   375 	  ++dfs;
   376 	  if (F.valid(typename MutableGraph::OutEdgeIt(dfs))) {
   377 	    if (dfs.isBNodeNewlyReached()) {
   378 	      typename MutableGraph::Node v=F.aNode(dfs);
   379 	      typename MutableGraph::Node w=F.bNode(dfs);
   380 	      pred.set(w, dfs);
   381 	      if (F.valid(pred.get(v))) {
   382 		free.set(w, std::min(free.get(v), residual_capacity.get(dfs)));
   383 	      } else {
   384 		free.set(w, residual_capacity.get(dfs)); 
   385 	      }
   386 	      if (w==tF) { 
   387 		__augment=true; 
   388 		_augment=true;
   389 		break; 
   390 	      }
   391 	      
   392 	    } else {
   393 	      F.erase(typename MutableGraph::OutEdgeIt(dfs));
   394 	    }
   395 	  } 
   396 	}
   397 
   398 	if (__augment) {
   399 	  typename MutableGraph::Node n=tF;
   400 	  Number augment_value=free.get(tF);
   401 	  while (F.valid(pred.get(n))) { 
   402 	    typename MutableGraph::Edge e=pred.get(n);
   403 	    res_graph.augment(original_edge.get(e), augment_value); 
   404 	    n=F.tail(e);
   405 	    if (residual_capacity.get(e)==augment_value) 
   406 	      F.erase(e); 
   407 	    else 
   408 	      residual_capacity.set(e, residual_capacity.get(e)-augment_value);
   409 	  }
   410 	}
   411 	
   412       }
   413             
   414       return _augment;
   415     }
   416     bool augmentOnBlockingFlow2() {
   417       bool _augment=false;
   418 
   419       //typedef ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> EAugGraph;
   420       typedef FilterGraphWrapper< ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> > EAugGraph;
   421       typedef typename EAugGraph::OutEdgeIt EAugOutEdgeIt;
   422       typedef typename EAugGraph::Edge EAugEdge;
   423 
   424       EAugGraph res_graph(*G, *flow, *capacity);
   425 
   426       //typedef typename EAugGraph::NodeMap<bool> ReachedMap;
   427       BfsIterator4< 
   428 	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>, 
   429 	typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt, 
   430 	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::NodeMap<bool> > bfs(res_graph);
   431       
   432       bfs.pushAndSetReached(s);
   433 
   434       typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::
   435 	NodeMap<int>& dist=res_graph.dist;
   436 
   437       while ( !bfs.finished() ) {
   438 	typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt e=bfs;
   439 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   440 	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
   441 	}
   442 	++bfs;	
   443       } //computing distances from s in the residual graph
   444 
   445       bool __augment=true;
   446 
   447       while (__augment) {
   448 
   449 	__augment=false;
   450 	//computing blocking flow with dfs
   451 	typedef typename EAugGraph::NodeMap<bool> BlockingReachedMap;
   452 	DfsIterator4< EAugGraph, EAugOutEdgeIt, BlockingReachedMap > 
   453 	  dfs(res_graph);
   454 	typename EAugGraph::NodeMap<EAugEdge> pred(res_graph); 
   455 	pred.set(s, EAugEdge(INVALID));
   456 	//invalid iterators for sources
   457 
   458 	typename EAugGraph::NodeMap<Number> free(res_graph);
   459 
   460 	dfs.pushAndSetReached(s);
   461 	while (!dfs.finished()) {
   462 	  ++dfs;
   463 	  if (res_graph.valid(EAugOutEdgeIt(dfs))) { 
   464 	    if (dfs.isBNodeNewlyReached()) {
   465 	  
   466 	      typename EAugGraph::Node v=res_graph.aNode(dfs);
   467 	      typename EAugGraph::Node w=res_graph.bNode(dfs);
   468 
   469 	      pred.set(w, EAugOutEdgeIt(dfs));
   470 	      if (res_graph.valid(pred.get(v))) {
   471 		free.set(w, std::min(free.get(v), res_graph.free(dfs)));
   472 	      } else {
   473 		free.set(w, res_graph.free(dfs)); 
   474 	      }
   475 	      
   476 	      if (w==t) { 
   477 		__augment=true; 
   478 		_augment=true;
   479 		break; 
   480 	      }
   481 	    } else {
   482 	      res_graph.erase(dfs);
   483 	    }
   484 	  } 
   485 
   486 	}
   487 
   488 	if (__augment) {
   489 	  typename EAugGraph::Node n=t;
   490 	  Number augment_value=free.get(t);
   491 	  while (res_graph.valid(pred.get(n))) { 
   492 	    EAugEdge e=pred.get(n);
   493 	    res_graph.augment(e, augment_value);
   494 	    n=res_graph.tail(e);
   495 	    if (res_graph.free(e)==0)
   496 	      res_graph.erase(e);
   497 	  }
   498 	}
   499       
   500       }
   501             
   502       return _augment;
   503     }
   504     void run() {
   505       //int num_of_augmentations=0;
   506       while (augmentOnShortestPath()) { 
   507 	//while (augmentOnBlockingFlow<MutableGraph>()) { 
   508 	//std::cout << ++num_of_augmentations << " ";
   509 	//std::cout<<std::endl;
   510       } 
   511     }
   512     template<typename MutableGraph> void run() {
   513       //int num_of_augmentations=0;
   514       //while (augmentOnShortestPath()) { 
   515 	while (augmentOnBlockingFlow<MutableGraph>()) { 
   516 	//std::cout << ++num_of_augmentations << " ";
   517 	//std::cout<<std::endl;
   518       } 
   519     }
   520     Number flowValue() { 
   521       Number a=0;
   522       OutEdgeIt e;
   523       for(G->/*getF*/first(e, s); G->valid(e); G->next(e)) {
   524 	a+=flow->get(e);
   525       }
   526       return a;
   527     }
   528   };
   529 
   530 
   531   template <typename Graph, typename Number, typename FlowMap, typename CapacityMap>
   532   class MaxMatching {
   533   public:
   534     typedef typename Graph::Node Node;
   535     typedef typename Graph::NodeIt NodeIt;
   536     typedef typename Graph::Edge Edge;
   537     typedef typename Graph::EdgeIt EdgeIt;
   538     typedef typename Graph::OutEdgeIt OutEdgeIt;
   539     typedef typename Graph::InEdgeIt InEdgeIt;
   540 
   541     typedef typename Graph::NodeMap<bool> SMap;
   542     typedef typename Graph::NodeMap<bool> TMap;
   543   private:
   544     const Graph* G;
   545     SMap* S;
   546     TMap* T;
   547     //Node s;
   548     //Node t;
   549     FlowMap* flow;
   550     const CapacityMap* capacity;
   551     typedef ResGraphWrapper<Graph, Number, FlowMap, CapacityMap > AugGraph;
   552     typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
   553     typedef typename AugGraph::Edge AugEdge;
   554     typename Graph::NodeMap<int> used; //0
   555 
   556   public:
   557     MaxMatching(const Graph& _G, SMap& _S, TMap& _T, FlowMap& _flow, const CapacityMap& _capacity) : 
   558       G(&_G), S(&_S), T(&_T), flow(&_flow), capacity(&_capacity), used(_G) { }
   559     bool augmentOnShortestPath() {
   560       AugGraph res_graph(*G, *flow, *capacity);
   561       bool _augment=false;
   562       
   563       typedef typename AugGraph::NodeMap<bool> ReachedMap;
   564       BfsIterator5< AugGraph, /*AugOutEdgeIt,*/ ReachedMap > res_bfs(res_graph);
   565       typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
   566       for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
   567 	if ((S->get(s)) && (used.get(s)<1) ) {
   568 	  //Number u=0;
   569 	  //for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
   570 	  //u+=flow->get(e);
   571 	  //if (u<1) {
   572 	    res_bfs.pushAndSetReached(s);
   573 	    pred.set(s, AugEdge(INVALID));
   574 	    //}
   575 	}
   576       }
   577       
   578       typename AugGraph::NodeMap<Number> free(res_graph);
   579 	
   580       Node n;
   581       //searching for augmenting path
   582       while ( !res_bfs.finished() ) { 
   583 	AugOutEdgeIt e=res_bfs;
   584 	if (res_graph.valid(e) && res_bfs.isBNodeNewlyReached()) {
   585 	  Node v=res_graph.tail(e);
   586 	  Node w=res_graph.head(e);
   587 	  pred.set(w, e);
   588 	  if (res_graph.valid(pred.get(v))) {
   589 	    free.set(w, std::min(free.get(v), res_graph.free(e)));
   590 	  } else {
   591 	    free.set(w, res_graph.free(e)); 
   592 	  }
   593 	  n=res_graph.head(e);
   594 	  if (T->get(n) && (used.get(n)<1) ) { 
   595 	    //Number u=0;
   596 	    //for(InEdgeIt f=G->template first<InEdgeIt>(n); G->valid(f); G->next(f))
   597 	    //u+=flow->get(f);
   598 	    //if (u<1) {
   599 	      _augment=true; 
   600 	      break; 
   601 	      //}
   602 	  }
   603 	}
   604 	
   605 	++res_bfs;
   606       } //end of searching augmenting path
   607 
   608       if (_augment) {
   609 	//Node n=t;
   610 	used.set(n, 1); //mind2 vegen jav
   611 	Number augment_value=free.get(n);
   612 	while (res_graph.valid(pred.get(n))) { 
   613 	  AugEdge e=pred.get(n);
   614 	  res_graph.augment(e, augment_value); 
   615 	  n=res_graph.tail(e);
   616 	}
   617 	used.set(n, 1); //mind2 vegen jav
   618       }
   619 
   620       return _augment;
   621     }
   622 
   623 //     template<typename MutableGraph> bool augmentOnBlockingFlow() {      
   624 //       bool _augment=false;
   625 
   626 //       AugGraph res_graph(*G, *flow, *capacity);
   627 
   628 //       typedef typename AugGraph::NodeMap<bool> ReachedMap;
   629 //       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > bfs(res_graph);
   630 
   631 
   632 
   633 
   634 
   635 //       //typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
   636 //       for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
   637 // 	if (S->get(s)) {
   638 // 	  Number u=0;
   639 // 	  for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
   640 // 	    u+=flow->get(e);
   641 // 	  if (u<1) {
   642 // 	    res_bfs.pushAndSetReached(s);
   643 // 	    //pred.set(s, AugEdge(INVALID));
   644 // 	  }
   645 // 	}
   646 //       }
   647 
   648 
   649 
   650 
   651 //       //bfs.pushAndSetReached(s);
   652 //       typename AugGraph::NodeMap<int> dist(res_graph); //filled up with 0's
   653 //       while ( !bfs.finished() ) { 
   654 // 	AugOutEdgeIt e=bfs;
   655 // 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   656 // 	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
   657 // 	}
   658 	
   659 // 	++bfs;
   660 //       } //computing distances from s in the residual graph
   661 
   662 //       MutableGraph F;
   663 //       typename AugGraph::NodeMap<typename MutableGraph::Node> 
   664 // 	res_graph_to_F(res_graph);
   665 //       for(typename AugGraph::NodeIt n=res_graph.template first<typename AugGraph::NodeIt>(); res_graph.valid(n); res_graph.next(n)) {
   666 // 	res_graph_to_F.set(n, F.addNode());
   667 //       }
   668       
   669 //       typename MutableGraph::Node sF=res_graph_to_F.get(s);
   670 //       typename MutableGraph::Node tF=res_graph_to_F.get(t);
   671 
   672 //       typename MutableGraph::EdgeMap<AugEdge> original_edge(F);
   673 //       typename MutableGraph::EdgeMap<Number> residual_capacity(F);
   674 
   675 //       //Making F to the graph containing the edges of the residual graph 
   676 //       //which are in some shortest paths
   677 //       for(typename AugGraph::EdgeIt e=res_graph.template first<typename AugGraph::EdgeIt>(); res_graph.valid(e); res_graph.next(e)) {
   678 // 	if (dist.get(res_graph.head(e))==dist.get(res_graph.tail(e))+1) {
   679 // 	  typename MutableGraph::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
   680 // 	  original_edge.update();
   681 // 	  original_edge.set(f, e);
   682 // 	  residual_capacity.update();
   683 // 	  residual_capacity.set(f, res_graph.free(e));
   684 // 	} 
   685 //       }
   686 
   687 //       bool __augment=true;
   688 
   689 //       while (__augment) {
   690 // 	__augment=false;
   691 // 	//computing blocking flow with dfs
   692 // 	typedef typename MutableGraph::NodeMap<bool> BlockingReachedMap;
   693 // 	DfsIterator4< MutableGraph, typename MutableGraph::OutEdgeIt, BlockingReachedMap > dfs(F);
   694 // 	typename MutableGraph::NodeMap<typename MutableGraph::Edge> pred(F);
   695 // 	pred.set(sF, typename MutableGraph::Edge(INVALID));
   696 // 	//invalid iterators for sources
   697 
   698 // 	typename MutableGraph::NodeMap<Number> free(F);
   699 
   700 // 	dfs.pushAndSetReached(sF);      
   701 // 	while (!dfs.finished()) {
   702 // 	  ++dfs;
   703 // 	  if (F.valid(typename MutableGraph::OutEdgeIt(dfs))) {
   704 // 	    if (dfs.isBNodeNewlyReached()) {
   705 // 	      typename MutableGraph::Node v=F.aNode(dfs);
   706 // 	      typename MutableGraph::Node w=F.bNode(dfs);
   707 // 	      pred.set(w, dfs);
   708 // 	      if (F.valid(pred.get(v))) {
   709 // 		free.set(w, std::min(free.get(v), residual_capacity.get(dfs)));
   710 // 	      } else {
   711 // 		free.set(w, residual_capacity.get(dfs)); 
   712 // 	      }
   713 // 	      if (w==tF) { 
   714 // 		__augment=true; 
   715 // 		_augment=true;
   716 // 		break; 
   717 // 	      }
   718 	      
   719 // 	    } else {
   720 // 	      F.erase(typename MutableGraph::OutEdgeIt(dfs));
   721 // 	    }
   722 // 	  } 
   723 // 	}
   724 
   725 // 	if (__augment) {
   726 // 	  typename MutableGraph::Node n=tF;
   727 // 	  Number augment_value=free.get(tF);
   728 // 	  while (F.valid(pred.get(n))) { 
   729 // 	    typename MutableGraph::Edge e=pred.get(n);
   730 // 	    res_graph.augment(original_edge.get(e), augment_value); 
   731 // 	    n=F.tail(e);
   732 // 	    if (residual_capacity.get(e)==augment_value) 
   733 // 	      F.erase(e); 
   734 // 	    else 
   735 // 	      residual_capacity.set(e, residual_capacity.get(e)-augment_value);
   736 // 	  }
   737 // 	}
   738 	
   739 //       }
   740             
   741 //       return _augment;
   742 //     }
   743     bool augmentOnBlockingFlow2() {
   744       bool _augment=false;
   745 
   746       //typedef ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> EAugGraph;
   747       typedef FilterGraphWrapper< ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> > EAugGraph;
   748       typedef typename EAugGraph::OutEdgeIt EAugOutEdgeIt;
   749       typedef typename EAugGraph::Edge EAugEdge;
   750 
   751       EAugGraph res_graph(*G, *flow, *capacity);
   752 
   753       //typedef typename EAugGraph::NodeMap<bool> ReachedMap;
   754       BfsIterator4< 
   755 	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>, 
   756 	typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt, 
   757 	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::NodeMap<bool> > bfs(res_graph);
   758 
   759 
   760       //typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
   761       for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
   762 	if (S->get(s)) {
   763 	  Number u=0;
   764 	  for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
   765 	    u+=flow->get(e);
   766 	  if (u<1) {
   767 	    bfs.pushAndSetReached(s);
   768 	    //pred.set(s, AugEdge(INVALID));
   769 	  }
   770 	}
   771       }
   772 
   773       
   774       //bfs.pushAndSetReached(s);
   775 
   776       typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::
   777 	NodeMap<int>& dist=res_graph.dist;
   778 
   779       while ( !bfs.finished() ) {
   780 	typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt e=bfs;
   781 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   782 	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
   783 	}
   784 	++bfs;	
   785       } //computing distances from s in the residual graph
   786 
   787       bool __augment=true;
   788 
   789       while (__augment) {
   790 
   791 	__augment=false;
   792 	//computing blocking flow with dfs
   793 	typedef typename EAugGraph::NodeMap<bool> BlockingReachedMap;
   794 	DfsIterator4< EAugGraph, EAugOutEdgeIt, BlockingReachedMap > 
   795 	  dfs(res_graph);
   796 	typename EAugGraph::NodeMap<EAugEdge> pred(res_graph, INVALID); 
   797 	//pred.set(s, EAugEdge(INVALID));
   798 	//invalid iterators for sources
   799 
   800 	typename EAugGraph::NodeMap<Number> free(res_graph);
   801 
   802 
   803 	//typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
   804       for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
   805 	if (S->get(s)) {
   806 	  Number u=0;
   807 	  for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
   808 	    u+=flow->get(e);
   809 	  if (u<1) {
   810 	    dfs.pushAndSetReached(s);
   811 	    //pred.set(s, AugEdge(INVALID));
   812 	  }
   813 	}
   814       }
   815 
   816 
   817 
   818       //dfs.pushAndSetReached(s);
   819       typename EAugGraph::Node n;
   820 	while (!dfs.finished()) {
   821 	  ++dfs;
   822 	  if (res_graph.valid(EAugOutEdgeIt(dfs))) { 
   823 	    if (dfs.isBNodeNewlyReached()) {
   824 	  
   825 	      typename EAugGraph::Node v=res_graph.aNode(dfs);
   826 	      typename EAugGraph::Node w=res_graph.bNode(dfs);
   827 
   828 	      pred.set(w, EAugOutEdgeIt(dfs));
   829 	      if (res_graph.valid(pred.get(v))) {
   830 		free.set(w, std::min(free.get(v), res_graph.free(dfs)));
   831 	      } else {
   832 		free.set(w, res_graph.free(dfs)); 
   833 	      }
   834 	     
   835 	      n=w;
   836 	      if (T->get(w)) {
   837 		Number u=0;
   838 		for(InEdgeIt f=G->template first<InEdgeIt>(n); G->valid(f); G->next(f))
   839 		  u+=flow->get(f);
   840 		if (u<1) {
   841 		  __augment=true; 
   842 		  _augment=true;
   843 		  break; 
   844 		}
   845 	      }
   846 	    } else {
   847 	      res_graph.erase(dfs);
   848 	    }
   849 	  } 
   850 
   851 	}
   852 
   853 	if (__augment) {
   854 	  // typename EAugGraph::Node n=t;
   855 	  Number augment_value=free.get(n);
   856 	  while (res_graph.valid(pred.get(n))) { 
   857 	    EAugEdge e=pred.get(n);
   858 	    res_graph.augment(e, augment_value);
   859 	    n=res_graph.tail(e);
   860 	    if (res_graph.free(e)==0)
   861 	      res_graph.erase(e);
   862 	  }
   863 	}
   864       
   865       }
   866             
   867       return _augment;
   868     }
   869     void run() {
   870       //int num_of_augmentations=0;
   871       while (augmentOnShortestPath()) { 
   872 	//while (augmentOnBlockingFlow<MutableGraph>()) { 
   873 	//std::cout << ++num_of_augmentations << " ";
   874 	//std::cout<<std::endl;
   875       } 
   876     }
   877 //     template<typename MutableGraph> void run() {
   878 //       //int num_of_augmentations=0;
   879 //       //while (augmentOnShortestPath()) { 
   880 // 	while (augmentOnBlockingFlow<MutableGraph>()) { 
   881 // 	//std::cout << ++num_of_augmentations << " ";
   882 // 	//std::cout<<std::endl;
   883 //       } 
   884 //     } 
   885     Number flowValue() { 
   886       Number a=0;
   887       EdgeIt e;
   888       for(G->/*getF*/first(e); G->valid(e); G->next(e)) {
   889 	a+=flow->get(e);
   890       }
   891       return a;
   892     }
   893   };
   894 
   895 
   896 
   897 
   898 
   899   
   900 //   template <typename Graph, typename Number, typename FlowMap, typename CapacityMap>
   901 //   class MaxFlow2 {
   902 //   public:
   903 //     typedef typename Graph::Node Node;
   904 //     typedef typename Graph::Edge Edge;
   905 //     typedef typename Graph::EdgeIt EdgeIt;
   906 //     typedef typename Graph::OutEdgeIt OutEdgeIt;
   907 //     typedef typename Graph::InEdgeIt InEdgeIt;
   908 //   private:
   909 //     const Graph& G;
   910 //     std::list<Node>& S;
   911 //     std::list<Node>& T;
   912 //     FlowMap& flow;
   913 //     const CapacityMap& capacity;
   914 //     typedef ResGraphWrapper<Graph, Number, FlowMap, CapacityMap > AugGraph;
   915 //     typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
   916 //     typedef typename AugGraph::Edge AugEdge;
   917 //     typename Graph::NodeMap<bool> SMap;
   918 //     typename Graph::NodeMap<bool> TMap;
   919 //   public:
   920 //     MaxFlow2(const Graph& _G, std::list<Node>& _S, std::list<Node>& _T, FlowMap& _flow, const CapacityMap& _capacity) : G(_G), S(_S), T(_T), flow(_flow), capacity(_capacity), SMap(_G), TMap(_G) { 
   921 //       for(typename std::list<Node>::const_iterator i=S.begin(); 
   922 // 	  i!=S.end(); ++i) { 
   923 // 	SMap.set(*i, true); 
   924 //       }
   925 //       for (typename std::list<Node>::const_iterator i=T.begin(); 
   926 // 	   i!=T.end(); ++i) { 
   927 // 	TMap.set(*i, true); 
   928 //       }
   929 //     }
   930 //     bool augment() {
   931 //       AugGraph res_graph(G, flow, capacity);
   932 //       bool _augment=false;
   933 //       Node reached_t_node;
   934       
   935 //       typedef typename AugGraph::NodeMap<bool> ReachedMap;
   936 //       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > res_bfs(res_graph);
   937 //       for(typename std::list<Node>::const_iterator i=S.begin(); 
   938 // 	  i!=S.end(); ++i) {
   939 // 	res_bfs.pushAndSetReached(*i);
   940 //       }
   941 //       //res_bfs.pushAndSetReached(s);
   942 	
   943 //       typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
   944 //       //filled up with invalid iterators
   945       
   946 //       typename AugGraph::NodeMap<Number> free(res_graph);
   947 	
   948 //       //searching for augmenting path
   949 //       while ( !res_bfs.finished() ) { 
   950 // 	AugOutEdgeIt e=/*AugOutEdgeIt*/(res_bfs);
   951 // 	if (e.valid() && res_bfs.isBNodeNewlyReached()) {
   952 // 	  Node v=res_graph.tail(e);
   953 // 	  Node w=res_graph.head(e);
   954 // 	  pred.set(w, e);
   955 // 	  if (pred.get(v).valid()) {
   956 // 	    free.set(w, std::min(free.get(v), e.free()));
   957 // 	  } else {
   958 // 	    free.set(w, e.free()); 
   959 // 	  }
   960 // 	  if (TMap.get(res_graph.head(e))) { 
   961 // 	    _augment=true; 
   962 // 	    reached_t_node=res_graph.head(e);
   963 // 	    break; 
   964 // 	  }
   965 // 	}
   966 	
   967 // 	++res_bfs;
   968 //       } //end of searching augmenting path
   969 
   970 //       if (_augment) {
   971 // 	Node n=reached_t_node;
   972 // 	Number augment_value=free.get(reached_t_node);
   973 // 	while (pred.get(n).valid()) { 
   974 // 	  AugEdge e=pred.get(n);
   975 // 	  e.augment(augment_value); 
   976 // 	  n=res_graph.tail(e);
   977 // 	}
   978 //       }
   979 
   980 //       return _augment;
   981 //     }
   982 //     void run() {
   983 //       while (augment()) { } 
   984 //     }
   985 //     Number flowValue() { 
   986 //       Number a=0;
   987 //       for(typename std::list<Node>::const_iterator i=S.begin(); 
   988 // 	  i!=S.end(); ++i) { 
   989 // 	for(OutEdgeIt e=G.template first<OutEdgeIt>(*i); e.valid(); ++e) {
   990 // 	  a+=flow.get(e);
   991 // 	}
   992 // 	for(InEdgeIt e=G.template first<InEdgeIt>(*i); e.valid(); ++e) {
   993 // 	  a-=flow.get(e);
   994 // 	}
   995 //       }
   996 //       return a;
   997 //     }
   998 //   };
   999 
  1000 
  1001 
  1002 } // namespace hugo
  1003 
  1004 #endif //EDMONDS_KARP_H