src/work/edmonds_karp.h
author marci
Wed, 17 Mar 2004 17:04:41 +0000
changeset 197 fff43d9c7110
parent 196 8a9b9360463e
child 198 5cec393baade
permissions -rw-r--r--
.
     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 
   555   public:
   556     MaxMatching(const Graph& _G, SMap& _S, TMap& _T, FlowMap& _flow, const CapacityMap& _capacity) : 
   557       G(&_G), S(&_S), T(&_T), flow(&_flow), capacity(&_capacity) { }
   558     bool augmentOnShortestPath() {
   559       AugGraph res_graph(*G, *flow, *capacity);
   560       bool _augment=false;
   561       
   562       typedef typename AugGraph::NodeMap<bool> ReachedMap;
   563       BfsIterator5< AugGraph, /*AugOutEdgeIt,*/ ReachedMap > res_bfs(res_graph);
   564       typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
   565       for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
   566 	if (S->get(s)) {
   567 	  Number u=0;
   568 	  for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
   569 	    u+=flow->get(e);
   570 	  if (u<1) {
   571 	    res_bfs.pushAndSetReached(s);
   572 	    pred.set(s, AugEdge(INVALID));
   573 	  }
   574 	}
   575       }
   576       
   577       typename AugGraph::NodeMap<Number> free(res_graph);
   578 	
   579       Node n;
   580       //searching for augmenting path
   581       while ( !res_bfs.finished() ) { 
   582 	AugOutEdgeIt e=res_bfs;
   583 	if (res_graph.valid(e) && res_bfs.isBNodeNewlyReached()) {
   584 	  Node v=res_graph.tail(e);
   585 	  Node w=res_graph.head(e);
   586 	  pred.set(w, e);
   587 	  if (res_graph.valid(pred.get(v))) {
   588 	    free.set(w, std::min(free.get(v), res_graph.free(e)));
   589 	  } else {
   590 	    free.set(w, res_graph.free(e)); 
   591 	  }
   592 	  n=res_graph.head(e);
   593 	  if (T->get(n)) { 
   594 	    Number u=0;
   595 	    for(InEdgeIt f=G->template first<InEdgeIt>(n); G->valid(f); G->next(f))
   596 	      u+=flow->get(f);
   597 	    if (u<1) {
   598 	      _augment=true; 
   599 	      break; 
   600 	    }
   601 	  }
   602 	}
   603 	
   604 	++res_bfs;
   605       } //end of searching augmenting path
   606 
   607       if (_augment) {
   608 	//Node n=t;
   609 	Number augment_value=free.get(n);
   610 	while (res_graph.valid(pred.get(n))) { 
   611 	  AugEdge e=pred.get(n);
   612 	  res_graph.augment(e, augment_value); 
   613 	  n=res_graph.tail(e);
   614 	}
   615       }
   616 
   617       return _augment;
   618     }
   619 
   620 //     template<typename MutableGraph> bool augmentOnBlockingFlow() {      
   621 //       bool _augment=false;
   622 
   623 //       AugGraph res_graph(*G, *flow, *capacity);
   624 
   625 //       typedef typename AugGraph::NodeMap<bool> ReachedMap;
   626 //       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > bfs(res_graph);
   627 
   628 
   629 
   630 
   631 
   632 //       //typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
   633 //       for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
   634 // 	if (S->get(s)) {
   635 // 	  Number u=0;
   636 // 	  for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
   637 // 	    u+=flow->get(e);
   638 // 	  if (u<1) {
   639 // 	    res_bfs.pushAndSetReached(s);
   640 // 	    //pred.set(s, AugEdge(INVALID));
   641 // 	  }
   642 // 	}
   643 //       }
   644 
   645 
   646 
   647 
   648 //       //bfs.pushAndSetReached(s);
   649 //       typename AugGraph::NodeMap<int> dist(res_graph); //filled up with 0's
   650 //       while ( !bfs.finished() ) { 
   651 // 	AugOutEdgeIt e=bfs;
   652 // 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   653 // 	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
   654 // 	}
   655 	
   656 // 	++bfs;
   657 //       } //computing distances from s in the residual graph
   658 
   659 //       MutableGraph F;
   660 //       typename AugGraph::NodeMap<typename MutableGraph::Node> 
   661 // 	res_graph_to_F(res_graph);
   662 //       for(typename AugGraph::NodeIt n=res_graph.template first<typename AugGraph::NodeIt>(); res_graph.valid(n); res_graph.next(n)) {
   663 // 	res_graph_to_F.set(n, F.addNode());
   664 //       }
   665       
   666 //       typename MutableGraph::Node sF=res_graph_to_F.get(s);
   667 //       typename MutableGraph::Node tF=res_graph_to_F.get(t);
   668 
   669 //       typename MutableGraph::EdgeMap<AugEdge> original_edge(F);
   670 //       typename MutableGraph::EdgeMap<Number> residual_capacity(F);
   671 
   672 //       //Making F to the graph containing the edges of the residual graph 
   673 //       //which are in some shortest paths
   674 //       for(typename AugGraph::EdgeIt e=res_graph.template first<typename AugGraph::EdgeIt>(); res_graph.valid(e); res_graph.next(e)) {
   675 // 	if (dist.get(res_graph.head(e))==dist.get(res_graph.tail(e))+1) {
   676 // 	  typename MutableGraph::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
   677 // 	  original_edge.update();
   678 // 	  original_edge.set(f, e);
   679 // 	  residual_capacity.update();
   680 // 	  residual_capacity.set(f, res_graph.free(e));
   681 // 	} 
   682 //       }
   683 
   684 //       bool __augment=true;
   685 
   686 //       while (__augment) {
   687 // 	__augment=false;
   688 // 	//computing blocking flow with dfs
   689 // 	typedef typename MutableGraph::NodeMap<bool> BlockingReachedMap;
   690 // 	DfsIterator4< MutableGraph, typename MutableGraph::OutEdgeIt, BlockingReachedMap > dfs(F);
   691 // 	typename MutableGraph::NodeMap<typename MutableGraph::Edge> pred(F);
   692 // 	pred.set(sF, typename MutableGraph::Edge(INVALID));
   693 // 	//invalid iterators for sources
   694 
   695 // 	typename MutableGraph::NodeMap<Number> free(F);
   696 
   697 // 	dfs.pushAndSetReached(sF);      
   698 // 	while (!dfs.finished()) {
   699 // 	  ++dfs;
   700 // 	  if (F.valid(typename MutableGraph::OutEdgeIt(dfs))) {
   701 // 	    if (dfs.isBNodeNewlyReached()) {
   702 // 	      typename MutableGraph::Node v=F.aNode(dfs);
   703 // 	      typename MutableGraph::Node w=F.bNode(dfs);
   704 // 	      pred.set(w, dfs);
   705 // 	      if (F.valid(pred.get(v))) {
   706 // 		free.set(w, std::min(free.get(v), residual_capacity.get(dfs)));
   707 // 	      } else {
   708 // 		free.set(w, residual_capacity.get(dfs)); 
   709 // 	      }
   710 // 	      if (w==tF) { 
   711 // 		__augment=true; 
   712 // 		_augment=true;
   713 // 		break; 
   714 // 	      }
   715 	      
   716 // 	    } else {
   717 // 	      F.erase(typename MutableGraph::OutEdgeIt(dfs));
   718 // 	    }
   719 // 	  } 
   720 // 	}
   721 
   722 // 	if (__augment) {
   723 // 	  typename MutableGraph::Node n=tF;
   724 // 	  Number augment_value=free.get(tF);
   725 // 	  while (F.valid(pred.get(n))) { 
   726 // 	    typename MutableGraph::Edge e=pred.get(n);
   727 // 	    res_graph.augment(original_edge.get(e), augment_value); 
   728 // 	    n=F.tail(e);
   729 // 	    if (residual_capacity.get(e)==augment_value) 
   730 // 	      F.erase(e); 
   731 // 	    else 
   732 // 	      residual_capacity.set(e, residual_capacity.get(e)-augment_value);
   733 // 	  }
   734 // 	}
   735 	
   736 //       }
   737             
   738 //       return _augment;
   739 //     }
   740     bool augmentOnBlockingFlow2() {
   741       bool _augment=false;
   742 
   743       //typedef ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> EAugGraph;
   744       typedef FilterGraphWrapper< ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> > EAugGraph;
   745       typedef typename EAugGraph::OutEdgeIt EAugOutEdgeIt;
   746       typedef typename EAugGraph::Edge EAugEdge;
   747 
   748       EAugGraph res_graph(*G, *flow, *capacity);
   749 
   750       //typedef typename EAugGraph::NodeMap<bool> ReachedMap;
   751       BfsIterator4< 
   752 	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>, 
   753 	typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt, 
   754 	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::NodeMap<bool> > bfs(res_graph);
   755 
   756 
   757       //typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
   758       for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
   759 	if (S->get(s)) {
   760 	  Number u=0;
   761 	  for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
   762 	    u+=flow->get(e);
   763 	  if (u<1) {
   764 	    bfs.pushAndSetReached(s);
   765 	    //pred.set(s, AugEdge(INVALID));
   766 	  }
   767 	}
   768       }
   769 
   770       
   771       //bfs.pushAndSetReached(s);
   772 
   773       typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::
   774 	NodeMap<int>& dist=res_graph.dist;
   775 
   776       while ( !bfs.finished() ) {
   777 	typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt e=bfs;
   778 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   779 	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
   780 	}
   781 	++bfs;	
   782       } //computing distances from s in the residual graph
   783 
   784       bool __augment=true;
   785 
   786       while (__augment) {
   787 
   788 	__augment=false;
   789 	//computing blocking flow with dfs
   790 	typedef typename EAugGraph::NodeMap<bool> BlockingReachedMap;
   791 	DfsIterator4< EAugGraph, EAugOutEdgeIt, BlockingReachedMap > 
   792 	  dfs(res_graph);
   793 	typename EAugGraph::NodeMap<EAugEdge> pred(res_graph, INVALID); 
   794 	//pred.set(s, EAugEdge(INVALID));
   795 	//invalid iterators for sources
   796 
   797 	typename EAugGraph::NodeMap<Number> free(res_graph);
   798 
   799 
   800 	//typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
   801       for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
   802 	if (S->get(s)) {
   803 	  Number u=0;
   804 	  for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
   805 	    u+=flow->get(e);
   806 	  if (u<1) {
   807 	    dfs.pushAndSetReached(s);
   808 	    //pred.set(s, AugEdge(INVALID));
   809 	  }
   810 	}
   811       }
   812 
   813 
   814 
   815       //dfs.pushAndSetReached(s);
   816       typename EAugGraph::Node n;
   817 	while (!dfs.finished()) {
   818 	  ++dfs;
   819 	  if (res_graph.valid(EAugOutEdgeIt(dfs))) { 
   820 	    if (dfs.isBNodeNewlyReached()) {
   821 	  
   822 	      typename EAugGraph::Node v=res_graph.aNode(dfs);
   823 	      typename EAugGraph::Node w=res_graph.bNode(dfs);
   824 
   825 	      pred.set(w, EAugOutEdgeIt(dfs));
   826 	      if (res_graph.valid(pred.get(v))) {
   827 		free.set(w, std::min(free.get(v), res_graph.free(dfs)));
   828 	      } else {
   829 		free.set(w, res_graph.free(dfs)); 
   830 	      }
   831 	     
   832 	      n=w;
   833 	      if (T->get(w)) {
   834 		Number u=0;
   835 		for(InEdgeIt f=G->template first<InEdgeIt>(n); G->valid(f); G->next(f))
   836 		  u+=flow->get(f);
   837 		if (u<1) {
   838 		  __augment=true; 
   839 		  _augment=true;
   840 		  break; 
   841 		}
   842 	      }
   843 	    } else {
   844 	      res_graph.erase(dfs);
   845 	    }
   846 	  } 
   847 
   848 	}
   849 
   850 	if (__augment) {
   851 	  // typename EAugGraph::Node n=t;
   852 	  Number augment_value=free.get(n);
   853 	  while (res_graph.valid(pred.get(n))) { 
   854 	    EAugEdge e=pred.get(n);
   855 	    res_graph.augment(e, augment_value);
   856 	    n=res_graph.tail(e);
   857 	    if (res_graph.free(e)==0)
   858 	      res_graph.erase(e);
   859 	  }
   860 	}
   861       
   862       }
   863             
   864       return _augment;
   865     }
   866     void run() {
   867       //int num_of_augmentations=0;
   868       while (augmentOnShortestPath()) { 
   869 	//while (augmentOnBlockingFlow<MutableGraph>()) { 
   870 	//std::cout << ++num_of_augmentations << " ";
   871 	//std::cout<<std::endl;
   872       } 
   873     }
   874 //     template<typename MutableGraph> void run() {
   875 //       //int num_of_augmentations=0;
   876 //       //while (augmentOnShortestPath()) { 
   877 // 	while (augmentOnBlockingFlow<MutableGraph>()) { 
   878 // 	//std::cout << ++num_of_augmentations << " ";
   879 // 	//std::cout<<std::endl;
   880 //       } 
   881 //     } 
   882     Number flowValue() { 
   883       Number a=0;
   884       EdgeIt e;
   885       for(G->/*getF*/first(e); G->valid(e); G->next(e)) {
   886 	a+=flow->get(e);
   887       }
   888       return a;
   889     }
   890   };
   891 
   892 
   893 
   894 
   895 
   896   
   897 //   template <typename Graph, typename Number, typename FlowMap, typename CapacityMap>
   898 //   class MaxFlow2 {
   899 //   public:
   900 //     typedef typename Graph::Node Node;
   901 //     typedef typename Graph::Edge Edge;
   902 //     typedef typename Graph::EdgeIt EdgeIt;
   903 //     typedef typename Graph::OutEdgeIt OutEdgeIt;
   904 //     typedef typename Graph::InEdgeIt InEdgeIt;
   905 //   private:
   906 //     const Graph& G;
   907 //     std::list<Node>& S;
   908 //     std::list<Node>& T;
   909 //     FlowMap& flow;
   910 //     const CapacityMap& capacity;
   911 //     typedef ResGraphWrapper<Graph, Number, FlowMap, CapacityMap > AugGraph;
   912 //     typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
   913 //     typedef typename AugGraph::Edge AugEdge;
   914 //     typename Graph::NodeMap<bool> SMap;
   915 //     typename Graph::NodeMap<bool> TMap;
   916 //   public:
   917 //     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) { 
   918 //       for(typename std::list<Node>::const_iterator i=S.begin(); 
   919 // 	  i!=S.end(); ++i) { 
   920 // 	SMap.set(*i, true); 
   921 //       }
   922 //       for (typename std::list<Node>::const_iterator i=T.begin(); 
   923 // 	   i!=T.end(); ++i) { 
   924 // 	TMap.set(*i, true); 
   925 //       }
   926 //     }
   927 //     bool augment() {
   928 //       AugGraph res_graph(G, flow, capacity);
   929 //       bool _augment=false;
   930 //       Node reached_t_node;
   931       
   932 //       typedef typename AugGraph::NodeMap<bool> ReachedMap;
   933 //       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > res_bfs(res_graph);
   934 //       for(typename std::list<Node>::const_iterator i=S.begin(); 
   935 // 	  i!=S.end(); ++i) {
   936 // 	res_bfs.pushAndSetReached(*i);
   937 //       }
   938 //       //res_bfs.pushAndSetReached(s);
   939 	
   940 //       typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
   941 //       //filled up with invalid iterators
   942       
   943 //       typename AugGraph::NodeMap<Number> free(res_graph);
   944 	
   945 //       //searching for augmenting path
   946 //       while ( !res_bfs.finished() ) { 
   947 // 	AugOutEdgeIt e=/*AugOutEdgeIt*/(res_bfs);
   948 // 	if (e.valid() && res_bfs.isBNodeNewlyReached()) {
   949 // 	  Node v=res_graph.tail(e);
   950 // 	  Node w=res_graph.head(e);
   951 // 	  pred.set(w, e);
   952 // 	  if (pred.get(v).valid()) {
   953 // 	    free.set(w, std::min(free.get(v), e.free()));
   954 // 	  } else {
   955 // 	    free.set(w, e.free()); 
   956 // 	  }
   957 // 	  if (TMap.get(res_graph.head(e))) { 
   958 // 	    _augment=true; 
   959 // 	    reached_t_node=res_graph.head(e);
   960 // 	    break; 
   961 // 	  }
   962 // 	}
   963 	
   964 // 	++res_bfs;
   965 //       } //end of searching augmenting path
   966 
   967 //       if (_augment) {
   968 // 	Node n=reached_t_node;
   969 // 	Number augment_value=free.get(reached_t_node);
   970 // 	while (pred.get(n).valid()) { 
   971 // 	  AugEdge e=pred.get(n);
   972 // 	  e.augment(augment_value); 
   973 // 	  n=res_graph.tail(e);
   974 // 	}
   975 //       }
   976 
   977 //       return _augment;
   978 //     }
   979 //     void run() {
   980 //       while (augment()) { } 
   981 //     }
   982 //     Number flowValue() { 
   983 //       Number a=0;
   984 //       for(typename std::list<Node>::const_iterator i=S.begin(); 
   985 // 	  i!=S.end(); ++i) { 
   986 // 	for(OutEdgeIt e=G.template first<OutEdgeIt>(*i); e.valid(); ++e) {
   987 // 	  a+=flow.get(e);
   988 // 	}
   989 // 	for(InEdgeIt e=G.template first<InEdgeIt>(*i); e.valid(); ++e) {
   990 // 	  a-=flow.get(e);
   991 // 	}
   992 //       }
   993 //       return a;
   994 //     }
   995 //   };
   996 
   997 
   998 
   999 } // namespace hugo
  1000 
  1001 #endif //EDMONDS_KARP_H