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