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