src/work/marci/experiment/edmonds_karp.h
changeset 281 3fefabfd00b7
child 921 818510fa3d99
equal deleted inserted replaced
-1:000000000000 0:b3fb35147b12
       
     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 
       
    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 GraphWrapper, typename Number, typename FlowMap, typename CapacityMap>
       
   251   class MaxFlow {
       
   252   protected:
       
   253     typedef GraphWrapper GW;
       
   254     typedef typename GW::Node Node;
       
   255     typedef typename GW::Edge Edge;
       
   256     typedef typename GW::EdgeIt EdgeIt;
       
   257     typedef typename GW::OutEdgeIt OutEdgeIt;
       
   258     typedef typename GW::InEdgeIt InEdgeIt;
       
   259     //const Graph* G;
       
   260     GW gw;
       
   261     Node s;
       
   262     Node t;
       
   263     FlowMap* flow;
       
   264     const CapacityMap* capacity;
       
   265     typedef ResGraphWrapper<GW, Number, FlowMap, CapacityMap > ResGW;
       
   266     typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
       
   267     typedef typename ResGW::Edge ResGWEdge;
       
   268   public:
       
   269 
       
   270     MaxFlow(const GW& _gw, Node _s, Node _t, FlowMap& _flow, const CapacityMap& _capacity) : 
       
   271       gw(_gw), s(_s), t(_t), flow(&_flow), capacity(&_capacity) { }
       
   272 
       
   273     bool augmentOnShortestPath() {
       
   274       ResGW res_graph(gw, *flow, *capacity);
       
   275       bool _augment=false;
       
   276       
       
   277       typedef typename ResGW::NodeMap<bool> ReachedMap;
       
   278       BfsIterator5< ResGW, ReachedMap > bfs(res_graph);
       
   279       bfs.pushAndSetReached(s);
       
   280 	
       
   281       typename ResGW::NodeMap<ResGWEdge> pred(res_graph); 
       
   282       pred.set(s, INVALID);
       
   283       
       
   284       typename ResGW::NodeMap<Number> free(res_graph);
       
   285 	
       
   286       //searching for augmenting path
       
   287       while ( !bfs.finished() ) { 
       
   288 	ResGWOutEdgeIt e=bfs;
       
   289 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
       
   290 	  Node v=res_graph.tail(e);
       
   291 	  Node w=res_graph.head(e);
       
   292 	  pred.set(w, e);
       
   293 	  if (res_graph.valid(pred.get(v))) {
       
   294 	    free.set(w, std::min(free.get(v), res_graph.resCap(e)));
       
   295 	  } else {
       
   296 	    free.set(w, res_graph.resCap(e)); 
       
   297 	  }
       
   298 	  if (res_graph.head(e)==t) { _augment=true; break; }
       
   299 	}
       
   300 	
       
   301 	++bfs;
       
   302       } //end of searching augmenting path
       
   303 
       
   304       if (_augment) {
       
   305 	Node n=t;
       
   306 	Number augment_value=free.get(t);
       
   307 	while (res_graph.valid(pred.get(n))) { 
       
   308 	  ResGWEdge e=pred.get(n);
       
   309 	  res_graph.augment(e, augment_value); 
       
   310 	  n=res_graph.tail(e);
       
   311 	}
       
   312       }
       
   313 
       
   314       return _augment;
       
   315     }
       
   316 
       
   317     template<typename MapGraphWrapper> 
       
   318     class DistanceMap {
       
   319     protected:
       
   320       MapGraphWrapper gw;
       
   321       typename MapGraphWrapper::NodeMap<int> dist; 
       
   322     public:
       
   323       DistanceMap(MapGraphWrapper& _gw) : gw(_gw), dist(_gw, _gw.nodeNum()) { }
       
   324       void set(const typename MapGraphWrapper::Node& n, int a) { dist[n]=a; }
       
   325       int get(const typename MapGraphWrapper::Node& n) const { return dist[n]; }
       
   326       bool get(const typename MapGraphWrapper::Edge& e) const { 
       
   327 	return (dist.get(gw.tail(e))<dist.get(gw.head(e))); 
       
   328       }
       
   329     };
       
   330 
       
   331     template<typename MutableGraph> bool augmentOnBlockingFlow() {      
       
   332       typedef MutableGraph MG;
       
   333       bool _augment=false;
       
   334 
       
   335       ResGW res_graph(gw, *flow, *capacity);
       
   336 
       
   337       typedef typename ResGW::NodeMap<bool> ReachedMap;
       
   338       BfsIterator5< ResGW, ReachedMap > bfs(res_graph);
       
   339 
       
   340       bfs.pushAndSetReached(s);
       
   341       //typename ResGW::NodeMap<int> dist(res_graph); //filled up with 0'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.get(res_graph.tail(e))+1);
       
   347 	}
       
   348 	++bfs;
       
   349       } //computing distances from s in the residual graph
       
   350 
       
   351       MG F;
       
   352       typedef SubGraphWrapper<ResGW, DistanceMap<ResGW> > FilterResGW;
       
   353       FilterResGW filter_res_graph(res_graph, dist);
       
   354       typename ResGW::NodeMap<typename MG::Node> res_graph_to_F(res_graph);
       
   355       {
       
   356 	typename ResGW::NodeIt n;
       
   357 	for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
       
   358 	  res_graph_to_F.set(n, F.addNode());
       
   359 	}
       
   360       }
       
   361 
       
   362       typename MG::Node sF=res_graph_to_F.get(s);
       
   363       typename MG::Node tF=res_graph_to_F.get(t);
       
   364       typename MG::EdgeMap<ResGWEdge> original_edge(F);
       
   365       typename MG::EdgeMap<Number> residual_capacity(F);
       
   366 
       
   367       //Making F to the graph containing the edges of the residual graph 
       
   368       //which are in some shortest paths
       
   369       {
       
   370 	typename FilterResGW::EdgeIt e;
       
   371 	for(filter_res_graph.first(e); filter_res_graph.valid(e); filter_res_graph.next(e)) {
       
   372 	  //if (dist.get(res_graph.head(e))==dist.get(res_graph.tail(e))+1) {
       
   373 	  typename MG::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
       
   374 	  original_edge.update();
       
   375 	  original_edge.set(f, e);
       
   376 	  residual_capacity.update();
       
   377 	  residual_capacity.set(f, res_graph.resCap(e));
       
   378 	  //} 
       
   379 	}
       
   380       }
       
   381 
       
   382       bool __augment=true;
       
   383 
       
   384       while (__augment) {
       
   385 	__augment=false;
       
   386 	//computing blocking flow with dfs
       
   387 	typedef typename TrivGraphWrapper<MG>::NodeMap<bool> BlockingReachedMap;
       
   388 	DfsIterator5< TrivGraphWrapper<MG>, BlockingReachedMap > dfs(F);
       
   389 	typename MG::NodeMap<typename MG::Edge> pred(F);
       
   390 	pred.set(sF, INVALID);
       
   391 	//invalid iterators for sources
       
   392 
       
   393 	typename MG::NodeMap<Number> free(F);
       
   394 
       
   395 	dfs.pushAndSetReached(sF);      
       
   396 	while (!dfs.finished()) {
       
   397 	  ++dfs;
       
   398 	  if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
       
   399 	    if (dfs.isBNodeNewlyReached()) {
       
   400 	      typename MG::Node v=F.aNode(dfs);
       
   401 	      typename MG::Node w=F.bNode(dfs);
       
   402 	      pred.set(w, dfs);
       
   403 	      if (F.valid(pred.get(v))) {
       
   404 		free.set(w, std::min(free.get(v), residual_capacity.get(dfs)));
       
   405 	      } else {
       
   406 		free.set(w, residual_capacity.get(dfs)); 
       
   407 	      }
       
   408 	      if (w==tF) { 
       
   409 		__augment=true; 
       
   410 		_augment=true;
       
   411 		break; 
       
   412 	      }
       
   413 	      
       
   414 	    } else {
       
   415 	      F.erase(/*typename MG::OutEdgeIt*/(dfs));
       
   416 	    }
       
   417 	  } 
       
   418 	}
       
   419 
       
   420 	if (__augment) {
       
   421 	  typename MG::Node n=tF;
       
   422 	  Number augment_value=free.get(tF);
       
   423 	  while (F.valid(pred.get(n))) { 
       
   424 	    typename MG::Edge e=pred.get(n);
       
   425 	    res_graph.augment(original_edge.get(e), augment_value); 
       
   426 	    n=F.tail(e);
       
   427 	    if (residual_capacity.get(e)==augment_value) 
       
   428 	      F.erase(e); 
       
   429 	    else 
       
   430 	      residual_capacity.set(e, residual_capacity.get(e)-augment_value);
       
   431 	  }
       
   432 	}
       
   433 	
       
   434       }
       
   435             
       
   436       return _augment;
       
   437     }
       
   438 
       
   439     template<typename MutableGraph> bool augmentOnBlockingFlow1() {      
       
   440       typedef MutableGraph MG;
       
   441       bool _augment=false;
       
   442 
       
   443       ResGW res_graph(gw, *flow, *capacity);
       
   444 
       
   445       //bfs for distances on the residual graph
       
   446       typedef typename ResGW::NodeMap<bool> ReachedMap;
       
   447       BfsIterator5< ResGW, ReachedMap > bfs(res_graph);
       
   448       bfs.pushAndSetReached(s);
       
   449       typename ResGW::NodeMap<int> dist(res_graph); //filled up with 0's
       
   450 
       
   451       //F will contain the physical copy of the residual graph
       
   452       //with the set of edges which are on shortest paths
       
   453       MG F;
       
   454       typename ResGW::NodeMap<typename MG::Node> res_graph_to_F(res_graph);
       
   455       {
       
   456 	typename ResGW::NodeIt n;
       
   457 	for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
       
   458 	  res_graph_to_F.set(n, F.addNode());
       
   459 	}
       
   460       }
       
   461 
       
   462       typename MG::Node sF=res_graph_to_F.get(s);
       
   463       typename MG::Node tF=res_graph_to_F.get(t);
       
   464       typename MG::EdgeMap<ResGWEdge> original_edge(F);
       
   465       typename MG::EdgeMap<Number> residual_capacity(F);
       
   466 
       
   467       while ( !bfs.finished() ) { 
       
   468 	ResGWOutEdgeIt e=bfs;
       
   469 	if (res_graph.valid(e)) {
       
   470 	  if (bfs.isBNodeNewlyReached()) {
       
   471 	    dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
       
   472 	    typename MG::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
       
   473 	    original_edge.update();
       
   474 	    original_edge.set(f, e);
       
   475 	    residual_capacity.update();
       
   476 	    residual_capacity.set(f, res_graph.resCap(e));
       
   477 	  } else {
       
   478 	    if (dist.get(res_graph.head(e))==(dist.get(res_graph.tail(e))+1)) {
       
   479 	      typename MG::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
       
   480 	      original_edge.update();
       
   481 	      original_edge.set(f, e);
       
   482 	      residual_capacity.update();
       
   483 	      residual_capacity.set(f, res_graph.resCap(e));
       
   484 	    }
       
   485 	  }
       
   486 	}
       
   487 	++bfs;
       
   488       } //computing distances from s in the residual graph
       
   489 
       
   490       bool __augment=true;
       
   491 
       
   492       while (__augment) {
       
   493 	__augment=false;
       
   494 	//computing blocking flow with dfs
       
   495 	typedef typename TrivGraphWrapper<MG>::NodeMap<bool> BlockingReachedMap;
       
   496 	DfsIterator5< TrivGraphWrapper<MG>, BlockingReachedMap > dfs(F);
       
   497 	typename MG::NodeMap<typename MG::Edge> pred(F);
       
   498 	pred.set(sF, INVALID);
       
   499 	//invalid iterators for sources
       
   500 
       
   501 	typename MG::NodeMap<Number> free(F);
       
   502 
       
   503 	dfs.pushAndSetReached(sF);      
       
   504 	while (!dfs.finished()) {
       
   505 	  ++dfs;
       
   506 	  if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
       
   507 	    if (dfs.isBNodeNewlyReached()) {
       
   508 	      typename MG::Node v=F.aNode(dfs);
       
   509 	      typename MG::Node w=F.bNode(dfs);
       
   510 	      pred.set(w, dfs);
       
   511 	      if (F.valid(pred.get(v))) {
       
   512 		free.set(w, std::min(free.get(v), residual_capacity.get(dfs)));
       
   513 	      } else {
       
   514 		free.set(w, residual_capacity.get(dfs)); 
       
   515 	      }
       
   516 	      if (w==tF) { 
       
   517 		__augment=true; 
       
   518 		_augment=true;
       
   519 		break; 
       
   520 	      }
       
   521 	      
       
   522 	    } else {
       
   523 	      F.erase(/*typename MG::OutEdgeIt*/(dfs));
       
   524 	    }
       
   525 	  } 
       
   526 	}
       
   527 
       
   528 	if (__augment) {
       
   529 	  typename MG::Node n=tF;
       
   530 	  Number augment_value=free.get(tF);
       
   531 	  while (F.valid(pred.get(n))) { 
       
   532 	    typename MG::Edge e=pred.get(n);
       
   533 	    res_graph.augment(original_edge.get(e), augment_value); 
       
   534 	    n=F.tail(e);
       
   535 	    if (residual_capacity.get(e)==augment_value) 
       
   536 	      F.erase(e); 
       
   537 	    else 
       
   538 	      residual_capacity.set(e, residual_capacity.get(e)-augment_value);
       
   539 	  }
       
   540 	}
       
   541 	
       
   542       }
       
   543             
       
   544       return _augment;
       
   545     }
       
   546 
       
   547     bool augmentOnBlockingFlow2() {
       
   548       bool _augment=false;
       
   549 
       
   550       ResGW res_graph(gw, *flow, *capacity);
       
   551 
       
   552       typedef typename ResGW::NodeMap<bool> ReachedMap;
       
   553       BfsIterator5< ResGW, ReachedMap > bfs(res_graph);
       
   554 
       
   555       bfs.pushAndSetReached(s);
       
   556       DistanceMap<ResGW> dist(res_graph);
       
   557       while ( !bfs.finished() ) { 
       
   558  	ResGWOutEdgeIt e=bfs;
       
   559  	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
       
   560  	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
       
   561  	}
       
   562 	++bfs;
       
   563       } //computing distances from s in the residual graph
       
   564 
       
   565       //Subgraph containing the edges on some shortest paths
       
   566       typedef SubGraphWrapper<ResGW, DistanceMap<ResGW> > FilterResGW;
       
   567       FilterResGW filter_res_graph(res_graph, dist);
       
   568 
       
   569       //Subgraph, which is able to delete edges which are already 
       
   570       //met by the dfs
       
   571       typename FilterResGW::NodeMap<typename FilterResGW::OutEdgeIt> 
       
   572  	first_out_edges(filter_res_graph);
       
   573       typename FilterResGW::NodeIt v;
       
   574       for(filter_res_graph.first(v); filter_res_graph.valid(v); 
       
   575  	  filter_res_graph.next(v)) 
       
   576       {
       
   577  	typename FilterResGW::OutEdgeIt e;
       
   578  	filter_res_graph.first(e, v);
       
   579  	first_out_edges.set(v, e);
       
   580       }
       
   581       typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
       
   582 	NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
       
   583       ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
       
   584 
       
   585       bool __augment=true;
       
   586 
       
   587       while (__augment) {
       
   588 
       
   589  	__augment=false;
       
   590  	//computing blocking flow with dfs
       
   591 	typedef typename ErasingResGW::NodeMap<bool> BlockingReachedMap;
       
   592  	DfsIterator5< ErasingResGW, BlockingReachedMap > 
       
   593  	  dfs(erasing_res_graph);
       
   594  	typename ErasingResGW::NodeMap<typename ErasingResGW::OutEdgeIt> 
       
   595  	  pred(erasing_res_graph); 
       
   596  	pred.set(s, INVALID);
       
   597  	//invalid iterators for sources
       
   598 
       
   599  	typename ErasingResGW::NodeMap<Number> free(erasing_res_graph);
       
   600 
       
   601  	dfs.pushAndSetReached(s);
       
   602  	while (!dfs.finished()) {
       
   603  	  ++dfs;
       
   604  	  if (erasing_res_graph.valid(
       
   605  		/*typename ErasingResGW::OutEdgeIt*/(dfs))) 
       
   606  	  { 
       
   607  	    if (dfs.isBNodeNewlyReached()) {
       
   608 	  
       
   609  	      typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
       
   610  	      typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
       
   611 
       
   612  	      pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
       
   613  	      if (erasing_res_graph.valid(pred.get(v))) {
       
   614  		free.set(w, std::min(free.get(v), res_graph.resCap(dfs)));
       
   615  	      } else {
       
   616  		free.set(w, res_graph.resCap(dfs)); 
       
   617  	      }
       
   618 	      
       
   619  	      if (w==t) { 
       
   620  		__augment=true; 
       
   621  		_augment=true;
       
   622  		break; 
       
   623  	      }
       
   624 	    } else {
       
   625 	      erasing_res_graph.erase(dfs);
       
   626 	    }
       
   627 	  }
       
   628 	}	
       
   629 
       
   630  	if (__augment) {
       
   631  	  typename ErasingResGW::Node n=t;
       
   632  	  Number augment_value=free.get(n);
       
   633  	  while (erasing_res_graph.valid(pred.get(n))) { 
       
   634  	    typename ErasingResGW::OutEdgeIt e=pred.get(n);
       
   635  	    res_graph.augment(e, augment_value);
       
   636  	    n=erasing_res_graph.tail(e);
       
   637  	    if (res_graph.resCap(e)==0)
       
   638  	      erasing_res_graph.erase(e);
       
   639  	  }
       
   640  	}
       
   641       
       
   642       } //while (__augment) 
       
   643             
       
   644       return _augment;
       
   645     }
       
   646 
       
   647 //     bool augmentOnBlockingFlow2() {
       
   648 //       bool _augment=false;
       
   649 
       
   650 //       //typedef ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> EAugGraph;
       
   651 //       typedef FilterGraphWrapper< ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> > EAugGraph;
       
   652 //       typedef typename EAugGraph::OutEdgeIt EAugOutEdgeIt;
       
   653 //       typedef typename EAugGraph::Edge EAugEdge;
       
   654 
       
   655 //       EAugGraph res_graph(*G, *flow, *capacity);
       
   656 
       
   657 //       //typedef typename EAugGraph::NodeMap<bool> ReachedMap;
       
   658 //       BfsIterator5< 
       
   659 // 	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>, 
       
   660 // 	/*typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt,*/ 
       
   661 // 	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::NodeMap<bool> > bfs(res_graph);
       
   662       
       
   663 //       bfs.pushAndSetReached(s);
       
   664 
       
   665 //       typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::
       
   666 // 	NodeMap<int>& dist=res_graph.dist;
       
   667 
       
   668 //       while ( !bfs.finished() ) {
       
   669 // 	typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt e=bfs;
       
   670 // 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
       
   671 // 	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
       
   672 // 	}
       
   673 // 	++bfs;	
       
   674 //       } //computing distances from s in the residual graph
       
   675 
       
   676 //       bool __augment=true;
       
   677 
       
   678 //       while (__augment) {
       
   679 
       
   680 // 	__augment=false;
       
   681 // 	//computing blocking flow with dfs
       
   682 // 	typedef typename EAugGraph::NodeMap<bool> BlockingReachedMap;
       
   683 // 	DfsIterator5< EAugGraph/*, EAugOutEdgeIt*/, BlockingReachedMap > 
       
   684 // 	  dfs(res_graph);
       
   685 // 	typename EAugGraph::NodeMap<EAugEdge> pred(res_graph); 
       
   686 // 	pred.set(s, EAugEdge(INVALID));
       
   687 // 	//invalid iterators for sources
       
   688 
       
   689 // 	typename EAugGraph::NodeMap<Number> free(res_graph);
       
   690 
       
   691 // 	dfs.pushAndSetReached(s);
       
   692 // 	while (!dfs.finished()) {
       
   693 // 	  ++dfs;
       
   694 // 	  if (res_graph.valid(EAugOutEdgeIt(dfs))) { 
       
   695 // 	    if (dfs.isBNodeNewlyReached()) {
       
   696 	  
       
   697 // 	      typename EAugGraph::Node v=res_graph.aNode(dfs);
       
   698 // 	      typename EAugGraph::Node w=res_graph.bNode(dfs);
       
   699 
       
   700 // 	      pred.set(w, EAugOutEdgeIt(dfs));
       
   701 // 	      if (res_graph.valid(pred.get(v))) {
       
   702 // 		free.set(w, std::min(free.get(v), res_graph.free(dfs)));
       
   703 // 	      } else {
       
   704 // 		free.set(w, res_graph.free(dfs)); 
       
   705 // 	      }
       
   706 	      
       
   707 // 	      if (w==t) { 
       
   708 // 		__augment=true; 
       
   709 // 		_augment=true;
       
   710 // 		break; 
       
   711 // 	      }
       
   712 // 	    } else {
       
   713 // 	      res_graph.erase(dfs);
       
   714 // 	    }
       
   715 // 	  } 
       
   716 
       
   717 // 	}
       
   718 
       
   719 // 	if (__augment) {
       
   720 // 	  typename EAugGraph::Node n=t;
       
   721 // 	  Number augment_value=free.get(t);
       
   722 // 	  while (res_graph.valid(pred.get(n))) { 
       
   723 // 	    EAugEdge e=pred.get(n);
       
   724 // 	    res_graph.augment(e, augment_value);
       
   725 // 	    n=res_graph.tail(e);
       
   726 // 	    if (res_graph.free(e)==0)
       
   727 // 	      res_graph.erase(e);
       
   728 // 	  }
       
   729 // 	}
       
   730       
       
   731 //       }
       
   732             
       
   733 //       return _augment;
       
   734 //     }
       
   735 
       
   736     void run() {
       
   737       //int num_of_augmentations=0;
       
   738       while (augmentOnShortestPath()) { 
       
   739 	//while (augmentOnBlockingFlow<MutableGraph>()) { 
       
   740 	//std::cout << ++num_of_augmentations << " ";
       
   741 	//std::cout<<std::endl;
       
   742       } 
       
   743     }
       
   744 
       
   745     template<typename MutableGraph> void run() {
       
   746       //int num_of_augmentations=0;
       
   747       //while (augmentOnShortestPath()) { 
       
   748 	while (augmentOnBlockingFlow<MutableGraph>()) { 
       
   749 	//std::cout << ++num_of_augmentations << " ";
       
   750 	//std::cout<<std::endl;
       
   751       } 
       
   752     }
       
   753 
       
   754     Number flowValue() { 
       
   755       Number a=0;
       
   756       OutEdgeIt e;
       
   757       for(gw.first(e, s); gw.valid(e); gw.next(e)) {
       
   758 	a+=flow->get(e);
       
   759       }
       
   760       return a;
       
   761     }
       
   762 
       
   763   };
       
   764 
       
   765 
       
   766 //   template <typename Graph, typename Number, typename FlowMap, typename CapacityMap>
       
   767 //   class MaxMatching {
       
   768 //   public:
       
   769 //     typedef typename Graph::Node Node;
       
   770 //     typedef typename Graph::NodeIt NodeIt;
       
   771 //     typedef typename Graph::Edge Edge;
       
   772 //     typedef typename Graph::EdgeIt EdgeIt;
       
   773 //     typedef typename Graph::OutEdgeIt OutEdgeIt;
       
   774 //     typedef typename Graph::InEdgeIt InEdgeIt;
       
   775 
       
   776 //     typedef typename Graph::NodeMap<bool> SMap;
       
   777 //     typedef typename Graph::NodeMap<bool> TMap;
       
   778 //   private:
       
   779 //     const Graph* G;
       
   780 //     SMap* S;
       
   781 //     TMap* T;
       
   782 //     //Node s;
       
   783 //     //Node t;
       
   784 //     FlowMap* flow;
       
   785 //     const CapacityMap* capacity;
       
   786 //     typedef ResGraphWrapper<Graph, Number, FlowMap, CapacityMap > AugGraph;
       
   787 //     typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
       
   788 //     typedef typename AugGraph::Edge AugEdge;
       
   789 //     typename Graph::NodeMap<int> used; //0
       
   790 
       
   791 //   public:
       
   792 //     MaxMatching(const Graph& _G, SMap& _S, TMap& _T, FlowMap& _flow, const CapacityMap& _capacity) : 
       
   793 //       G(&_G), S(&_S), T(&_T), flow(&_flow), capacity(&_capacity), used(_G) { }
       
   794 //     bool augmentOnShortestPath() {
       
   795 //       AugGraph res_graph(*G, *flow, *capacity);
       
   796 //       bool _augment=false;
       
   797       
       
   798 //       typedef typename AugGraph::NodeMap<bool> ReachedMap;
       
   799 //       BfsIterator5< AugGraph, /*AugOutEdgeIt,*/ ReachedMap > bfs(res_graph);
       
   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)) && (used.get(s)<1) ) {
       
   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 // 	    bfs.pushAndSetReached(s);
       
   808 // 	    pred.set(s, AugEdge(INVALID));
       
   809 // 	    //}
       
   810 // 	}
       
   811 //       }
       
   812       
       
   813 //       typename AugGraph::NodeMap<Number> free(res_graph);
       
   814 	
       
   815 //       Node n;
       
   816 //       //searching for augmenting path
       
   817 //       while ( !bfs.finished() ) { 
       
   818 // 	AugOutEdgeIt e=bfs;
       
   819 // 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
       
   820 // 	  Node v=res_graph.tail(e);
       
   821 // 	  Node w=res_graph.head(e);
       
   822 // 	  pred.set(w, e);
       
   823 // 	  if (res_graph.valid(pred.get(v))) {
       
   824 // 	    free.set(w, std::min(free.get(v), res_graph.free(e)));
       
   825 // 	  } else {
       
   826 // 	    free.set(w, res_graph.free(e)); 
       
   827 // 	  }
       
   828 // 	  n=res_graph.head(e);
       
   829 // 	  if (T->get(n) && (used.get(n)<1) ) { 
       
   830 // 	    //Number u=0;
       
   831 // 	    //for(InEdgeIt f=G->template first<InEdgeIt>(n); G->valid(f); G->next(f))
       
   832 // 	    //u+=flow->get(f);
       
   833 // 	    //if (u<1) {
       
   834 // 	      _augment=true; 
       
   835 // 	      break; 
       
   836 // 	      //}
       
   837 // 	  }
       
   838 // 	}
       
   839 	
       
   840 // 	++bfs;
       
   841 //       } //end of searching augmenting path
       
   842 
       
   843 //       if (_augment) {
       
   844 // 	//Node n=t;
       
   845 // 	used.set(n, 1); //mind2 vegen jav
       
   846 // 	Number augment_value=free.get(n);
       
   847 // 	while (res_graph.valid(pred.get(n))) { 
       
   848 // 	  AugEdge e=pred.get(n);
       
   849 // 	  res_graph.augment(e, augment_value); 
       
   850 // 	  n=res_graph.tail(e);
       
   851 // 	}
       
   852 // 	used.set(n, 1); //mind2 vegen jav
       
   853 //       }
       
   854 
       
   855 //       return _augment;
       
   856 //     }
       
   857 
       
   858 // //     template<typename MutableGraph> bool augmentOnBlockingFlow() {      
       
   859 // //       bool _augment=false;
       
   860 
       
   861 // //       AugGraph res_graph(*G, *flow, *capacity);
       
   862 
       
   863 // //       typedef typename AugGraph::NodeMap<bool> ReachedMap;
       
   864 // //       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > bfs(res_graph);
       
   865 
       
   866 
       
   867 
       
   868 
       
   869 
       
   870 // //       //typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
       
   871 // //       for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
       
   872 // // 	if (S->get(s)) {
       
   873 // // 	  Number u=0;
       
   874 // // 	  for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
       
   875 // // 	    u+=flow->get(e);
       
   876 // // 	  if (u<1) {
       
   877 // // 	    bfs.pushAndSetReached(s);
       
   878 // // 	    //pred.set(s, AugEdge(INVALID));
       
   879 // // 	  }
       
   880 // // 	}
       
   881 // //       }
       
   882 
       
   883 
       
   884 
       
   885 
       
   886 // //       //bfs.pushAndSetReached(s);
       
   887 // //       typename AugGraph::NodeMap<int> dist(res_graph); //filled up with 0's
       
   888 // //       while ( !bfs.finished() ) { 
       
   889 // // 	AugOutEdgeIt e=bfs;
       
   890 // // 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
       
   891 // // 	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
       
   892 // // 	}
       
   893 	
       
   894 // // 	++bfs;
       
   895 // //       } //computing distances from s in the residual graph
       
   896 
       
   897 // //       MutableGraph F;
       
   898 // //       typename AugGraph::NodeMap<typename MutableGraph::Node> 
       
   899 // // 	res_graph_to_F(res_graph);
       
   900 // //       for(typename AugGraph::NodeIt n=res_graph.template first<typename AugGraph::NodeIt>(); res_graph.valid(n); res_graph.next(n)) {
       
   901 // // 	res_graph_to_F.set(n, F.addNode());
       
   902 // //       }
       
   903       
       
   904 // //       typename MutableGraph::Node sF=res_graph_to_F.get(s);
       
   905 // //       typename MutableGraph::Node tF=res_graph_to_F.get(t);
       
   906 
       
   907 // //       typename MutableGraph::EdgeMap<AugEdge> original_edge(F);
       
   908 // //       typename MutableGraph::EdgeMap<Number> residual_capacity(F);
       
   909 
       
   910 // //       //Making F to the graph containing the edges of the residual graph 
       
   911 // //       //which are in some shortest paths
       
   912 // //       for(typename AugGraph::EdgeIt e=res_graph.template first<typename AugGraph::EdgeIt>(); res_graph.valid(e); res_graph.next(e)) {
       
   913 // // 	if (dist.get(res_graph.head(e))==dist.get(res_graph.tail(e))+1) {
       
   914 // // 	  typename MutableGraph::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
       
   915 // // 	  original_edge.update();
       
   916 // // 	  original_edge.set(f, e);
       
   917 // // 	  residual_capacity.update();
       
   918 // // 	  residual_capacity.set(f, res_graph.free(e));
       
   919 // // 	} 
       
   920 // //       }
       
   921 
       
   922 // //       bool __augment=true;
       
   923 
       
   924 // //       while (__augment) {
       
   925 // // 	__augment=false;
       
   926 // // 	//computing blocking flow with dfs
       
   927 // // 	typedef typename MutableGraph::NodeMap<bool> BlockingReachedMap;
       
   928 // // 	DfsIterator4< MutableGraph, typename MutableGraph::OutEdgeIt, BlockingReachedMap > dfs(F);
       
   929 // // 	typename MutableGraph::NodeMap<typename MutableGraph::Edge> pred(F);
       
   930 // // 	pred.set(sF, typename MutableGraph::Edge(INVALID));
       
   931 // // 	//invalid iterators for sources
       
   932 
       
   933 // // 	typename MutableGraph::NodeMap<Number> free(F);
       
   934 
       
   935 // // 	dfs.pushAndSetReached(sF);      
       
   936 // // 	while (!dfs.finished()) {
       
   937 // // 	  ++dfs;
       
   938 // // 	  if (F.valid(typename MutableGraph::OutEdgeIt(dfs))) {
       
   939 // // 	    if (dfs.isBNodeNewlyReached()) {
       
   940 // // 	      typename MutableGraph::Node v=F.aNode(dfs);
       
   941 // // 	      typename MutableGraph::Node w=F.bNode(dfs);
       
   942 // // 	      pred.set(w, dfs);
       
   943 // // 	      if (F.valid(pred.get(v))) {
       
   944 // // 		free.set(w, std::min(free.get(v), residual_capacity.get(dfs)));
       
   945 // // 	      } else {
       
   946 // // 		free.set(w, residual_capacity.get(dfs)); 
       
   947 // // 	      }
       
   948 // // 	      if (w==tF) { 
       
   949 // // 		__augment=true; 
       
   950 // // 		_augment=true;
       
   951 // // 		break; 
       
   952 // // 	      }
       
   953 	      
       
   954 // // 	    } else {
       
   955 // // 	      F.erase(typename MutableGraph::OutEdgeIt(dfs));
       
   956 // // 	    }
       
   957 // // 	  } 
       
   958 // // 	}
       
   959 
       
   960 // // 	if (__augment) {
       
   961 // // 	  typename MutableGraph::Node n=tF;
       
   962 // // 	  Number augment_value=free.get(tF);
       
   963 // // 	  while (F.valid(pred.get(n))) { 
       
   964 // // 	    typename MutableGraph::Edge e=pred.get(n);
       
   965 // // 	    res_graph.augment(original_edge.get(e), augment_value); 
       
   966 // // 	    n=F.tail(e);
       
   967 // // 	    if (residual_capacity.get(e)==augment_value) 
       
   968 // // 	      F.erase(e); 
       
   969 // // 	    else 
       
   970 // // 	      residual_capacity.set(e, residual_capacity.get(e)-augment_value);
       
   971 // // 	  }
       
   972 // // 	}
       
   973 	
       
   974 // //       }
       
   975             
       
   976 // //       return _augment;
       
   977 // //     }
       
   978 //     bool augmentOnBlockingFlow2() {
       
   979 //       bool _augment=false;
       
   980 
       
   981 //       //typedef ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> EAugGraph;
       
   982 //       typedef FilterGraphWrapper< ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> > EAugGraph;
       
   983 //       typedef typename EAugGraph::OutEdgeIt EAugOutEdgeIt;
       
   984 //       typedef typename EAugGraph::Edge EAugEdge;
       
   985 
       
   986 //       EAugGraph res_graph(*G, *flow, *capacity);
       
   987 
       
   988 //       //typedef typename EAugGraph::NodeMap<bool> ReachedMap;
       
   989 //       BfsIterator5< 
       
   990 // 	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>, 
       
   991 // 	/*typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt,*/ 
       
   992 // 	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::NodeMap<bool> > bfs(res_graph);
       
   993 
       
   994 
       
   995 //       //typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
       
   996 //       for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
       
   997 // 	if (S->get(s)) {
       
   998 // 	  Number u=0;
       
   999 // 	  for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
       
  1000 // 	    u+=flow->get(e);
       
  1001 // 	  if (u<1) {
       
  1002 // 	    bfs.pushAndSetReached(s);
       
  1003 // 	    //pred.set(s, AugEdge(INVALID));
       
  1004 // 	  }
       
  1005 // 	}
       
  1006 //       }
       
  1007 
       
  1008       
       
  1009 //       //bfs.pushAndSetReached(s);
       
  1010 
       
  1011 //       typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::
       
  1012 // 	NodeMap<int>& dist=res_graph.dist;
       
  1013 
       
  1014 //       while ( !bfs.finished() ) {
       
  1015 // 	typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt e=bfs;
       
  1016 // 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
       
  1017 // 	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
       
  1018 // 	}
       
  1019 // 	++bfs;	
       
  1020 //       } //computing distances from s in the residual graph
       
  1021 
       
  1022 //       bool __augment=true;
       
  1023 
       
  1024 //       while (__augment) {
       
  1025 
       
  1026 // 	__augment=false;
       
  1027 // 	//computing blocking flow with dfs
       
  1028 // 	typedef typename EAugGraph::NodeMap<bool> BlockingReachedMap;
       
  1029 // 	DfsIterator5< EAugGraph/*, EAugOutEdgeIt*/, BlockingReachedMap > 
       
  1030 // 	  dfs(res_graph);
       
  1031 // 	typename EAugGraph::NodeMap<EAugEdge> pred(res_graph, INVALID); 
       
  1032 // 	//pred.set(s, EAugEdge(INVALID));
       
  1033 // 	//invalid iterators for sources
       
  1034 
       
  1035 // 	typename EAugGraph::NodeMap<Number> free(res_graph);
       
  1036 
       
  1037 
       
  1038 // 	//typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
       
  1039 //       for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
       
  1040 // 	if (S->get(s)) {
       
  1041 // 	  Number u=0;
       
  1042 // 	  for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
       
  1043 // 	    u+=flow->get(e);
       
  1044 // 	  if (u<1) {
       
  1045 // 	    dfs.pushAndSetReached(s);
       
  1046 // 	    //pred.set(s, AugEdge(INVALID));
       
  1047 // 	  }
       
  1048 // 	}
       
  1049 //       }
       
  1050 
       
  1051 
       
  1052 
       
  1053 //       //dfs.pushAndSetReached(s);
       
  1054 //       typename EAugGraph::Node n;
       
  1055 // 	while (!dfs.finished()) {
       
  1056 // 	  ++dfs;
       
  1057 // 	  if (res_graph.valid(EAugOutEdgeIt(dfs))) { 
       
  1058 // 	    if (dfs.isBNodeNewlyReached()) {
       
  1059 	  
       
  1060 // 	      typename EAugGraph::Node v=res_graph.aNode(dfs);
       
  1061 // 	      typename EAugGraph::Node w=res_graph.bNode(dfs);
       
  1062 
       
  1063 // 	      pred.set(w, EAugOutEdgeIt(dfs));
       
  1064 // 	      if (res_graph.valid(pred.get(v))) {
       
  1065 // 		free.set(w, std::min(free.get(v), res_graph.free(dfs)));
       
  1066 // 	      } else {
       
  1067 // 		free.set(w, res_graph.free(dfs)); 
       
  1068 // 	      }
       
  1069 	     
       
  1070 // 	      n=w;
       
  1071 // 	      if (T->get(w)) {
       
  1072 // 		Number u=0;
       
  1073 // 		for(InEdgeIt f=G->template first<InEdgeIt>(n); G->valid(f); G->next(f))
       
  1074 // 		  u+=flow->get(f);
       
  1075 // 		if (u<1) {
       
  1076 // 		  __augment=true; 
       
  1077 // 		  _augment=true;
       
  1078 // 		  break; 
       
  1079 // 		}
       
  1080 // 	      }
       
  1081 // 	    } else {
       
  1082 // 	      res_graph.erase(dfs);
       
  1083 // 	    }
       
  1084 // 	  } 
       
  1085 
       
  1086 // 	}
       
  1087 
       
  1088 // 	if (__augment) {
       
  1089 // 	  // typename EAugGraph::Node n=t;
       
  1090 // 	  Number augment_value=free.get(n);
       
  1091 // 	  while (res_graph.valid(pred.get(n))) { 
       
  1092 // 	    EAugEdge e=pred.get(n);
       
  1093 // 	    res_graph.augment(e, augment_value);
       
  1094 // 	    n=res_graph.tail(e);
       
  1095 // 	    if (res_graph.free(e)==0)
       
  1096 // 	      res_graph.erase(e);
       
  1097 // 	  }
       
  1098 // 	}
       
  1099       
       
  1100 //       }
       
  1101             
       
  1102 //       return _augment;
       
  1103 //     }
       
  1104 //     void run() {
       
  1105 //       //int num_of_augmentations=0;
       
  1106 //       while (augmentOnShortestPath()) { 
       
  1107 // 	//while (augmentOnBlockingFlow<MutableGraph>()) { 
       
  1108 // 	//std::cout << ++num_of_augmentations << " ";
       
  1109 // 	//std::cout<<std::endl;
       
  1110 //       } 
       
  1111 //     }
       
  1112 // //     template<typename MutableGraph> void run() {
       
  1113 // //       //int num_of_augmentations=0;
       
  1114 // //       //while (augmentOnShortestPath()) { 
       
  1115 // // 	while (augmentOnBlockingFlow<MutableGraph>()) { 
       
  1116 // // 	//std::cout << ++num_of_augmentations << " ";
       
  1117 // // 	//std::cout<<std::endl;
       
  1118 // //       } 
       
  1119 // //     } 
       
  1120 //     Number flowValue() { 
       
  1121 //       Number a=0;
       
  1122 //       EdgeIt e;
       
  1123 //       for(G->/*getF*/first(e); G->valid(e); G->next(e)) {
       
  1124 // 	a+=flow->get(e);
       
  1125 //       }
       
  1126 //       return a;
       
  1127 //     }
       
  1128 //   };
       
  1129 
       
  1130 
       
  1131 
       
  1132 
       
  1133 
       
  1134   
       
  1135 // //   template <typename Graph, typename Number, typename FlowMap, typename CapacityMap>
       
  1136 // //   class MaxFlow2 {
       
  1137 // //   public:
       
  1138 // //     typedef typename Graph::Node Node;
       
  1139 // //     typedef typename Graph::Edge Edge;
       
  1140 // //     typedef typename Graph::EdgeIt EdgeIt;
       
  1141 // //     typedef typename Graph::OutEdgeIt OutEdgeIt;
       
  1142 // //     typedef typename Graph::InEdgeIt InEdgeIt;
       
  1143 // //   private:
       
  1144 // //     const Graph& G;
       
  1145 // //     std::list<Node>& S;
       
  1146 // //     std::list<Node>& T;
       
  1147 // //     FlowMap& flow;
       
  1148 // //     const CapacityMap& capacity;
       
  1149 // //     typedef ResGraphWrapper<Graph, Number, FlowMap, CapacityMap > AugGraph;
       
  1150 // //     typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
       
  1151 // //     typedef typename AugGraph::Edge AugEdge;
       
  1152 // //     typename Graph::NodeMap<bool> SMap;
       
  1153 // //     typename Graph::NodeMap<bool> TMap;
       
  1154 // //   public:
       
  1155 // //     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) { 
       
  1156 // //       for(typename std::list<Node>::const_iterator i=S.begin(); 
       
  1157 // // 	  i!=S.end(); ++i) { 
       
  1158 // // 	SMap.set(*i, true); 
       
  1159 // //       }
       
  1160 // //       for (typename std::list<Node>::const_iterator i=T.begin(); 
       
  1161 // // 	   i!=T.end(); ++i) { 
       
  1162 // // 	TMap.set(*i, true); 
       
  1163 // //       }
       
  1164 // //     }
       
  1165 // //     bool augment() {
       
  1166 // //       AugGraph res_graph(G, flow, capacity);
       
  1167 // //       bool _augment=false;
       
  1168 // //       Node reached_t_node;
       
  1169       
       
  1170 // //       typedef typename AugGraph::NodeMap<bool> ReachedMap;
       
  1171 // //       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > bfs(res_graph);
       
  1172 // //       for(typename std::list<Node>::const_iterator i=S.begin(); 
       
  1173 // // 	  i!=S.end(); ++i) {
       
  1174 // // 	bfs.pushAndSetReached(*i);
       
  1175 // //       }
       
  1176 // //       //bfs.pushAndSetReached(s);
       
  1177 	
       
  1178 // //       typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
       
  1179 // //       //filled up with invalid iterators
       
  1180       
       
  1181 // //       typename AugGraph::NodeMap<Number> free(res_graph);
       
  1182 	
       
  1183 // //       //searching for augmenting path
       
  1184 // //       while ( !bfs.finished() ) { 
       
  1185 // // 	AugOutEdgeIt e=/*AugOutEdgeIt*/(bfs);
       
  1186 // // 	if (e.valid() && bfs.isBNodeNewlyReached()) {
       
  1187 // // 	  Node v=res_graph.tail(e);
       
  1188 // // 	  Node w=res_graph.head(e);
       
  1189 // // 	  pred.set(w, e);
       
  1190 // // 	  if (pred.get(v).valid()) {
       
  1191 // // 	    free.set(w, std::min(free.get(v), e.free()));
       
  1192 // // 	  } else {
       
  1193 // // 	    free.set(w, e.free()); 
       
  1194 // // 	  }
       
  1195 // // 	  if (TMap.get(res_graph.head(e))) { 
       
  1196 // // 	    _augment=true; 
       
  1197 // // 	    reached_t_node=res_graph.head(e);
       
  1198 // // 	    break; 
       
  1199 // // 	  }
       
  1200 // // 	}
       
  1201 	
       
  1202 // // 	++bfs;
       
  1203 // //       } //end of searching augmenting path
       
  1204 
       
  1205 // //       if (_augment) {
       
  1206 // // 	Node n=reached_t_node;
       
  1207 // // 	Number augment_value=free.get(reached_t_node);
       
  1208 // // 	while (pred.get(n).valid()) { 
       
  1209 // // 	  AugEdge e=pred.get(n);
       
  1210 // // 	  e.augment(augment_value); 
       
  1211 // // 	  n=res_graph.tail(e);
       
  1212 // // 	}
       
  1213 // //       }
       
  1214 
       
  1215 // //       return _augment;
       
  1216 // //     }
       
  1217 // //     void run() {
       
  1218 // //       while (augment()) { } 
       
  1219 // //     }
       
  1220 // //     Number flowValue() { 
       
  1221 // //       Number a=0;
       
  1222 // //       for(typename std::list<Node>::const_iterator i=S.begin(); 
       
  1223 // // 	  i!=S.end(); ++i) { 
       
  1224 // // 	for(OutEdgeIt e=G.template first<OutEdgeIt>(*i); e.valid(); ++e) {
       
  1225 // // 	  a+=flow.get(e);
       
  1226 // // 	}
       
  1227 // // 	for(InEdgeIt e=G.template first<InEdgeIt>(*i); e.valid(); ++e) {
       
  1228 // // 	  a-=flow.get(e);
       
  1229 // // 	}
       
  1230 // //       }
       
  1231 // //       return a;
       
  1232 // //     }
       
  1233 // //   };
       
  1234 
       
  1235 
       
  1236 } // namespace hugo
       
  1237 
       
  1238 #endif //HUGO_EDMONDS_KARP_H