Bugfix, thanks to Janos.
1 #ifndef EDMONDS_KARP_HH
2 #define EDMONDS_KARP_HH
8 #include <bfs_iterator.hh>
9 //#include <time_measure.h>
13 template<typename Graph, typename Number, typename FlowMap, typename CapacityMap>
16 typedef typename Graph::NodeIt NodeIt;
17 typedef typename Graph::EachNodeIt EachNodeIt;
19 typedef typename Graph::SymEdgeIt OldSymEdgeIt;
22 const CapacityMap& capacity;
24 ResGraph(const Graph& _G, FlowMap& _flow,
25 const CapacityMap& _capacity) :
26 G(_G), flow(_flow), capacity(_capacity) { }
31 friend class OutEdgeIt;
34 friend class ResGraph<Graph, Number, FlowMap, CapacityMap>;
36 const ResGraph<Graph, Number, FlowMap, CapacityMap>* resG;
40 //EdgeIt(const EdgeIt& e) : resG(e.resG), sym(e.sym) { }
42 if (resG->G.aNode(sym)==resG->G.tail(sym)) {
43 return (resG->capacity.get(sym)-resG->flow.get(sym));
45 return (resG->flow.get(sym));
48 bool valid() const { return sym.valid(); }
49 void augment(Number a) const {
50 if (resG->G.aNode(sym)==resG->G.tail(sym)) {
51 resG->flow.set(sym, resG->flow.get(sym)+a);
54 resG->flow.set(sym, resG->flow.get(sym)-a);
60 class OutEdgeIt : public EdgeIt {
61 friend class ResGraph<Graph, Number, FlowMap, CapacityMap>;
64 //OutEdgeIt(const OutEdgeIt& e) { resG=e.resG; sym=e.sym; }
66 OutEdgeIt(const ResGraph<Graph, Number, FlowMap, CapacityMap>& _resG, NodeIt v) {
68 sym=resG->G.template first<OldSymEdgeIt>(v);
69 while( sym.valid() && !(free()>0) ) { ++sym; }
72 OutEdgeIt& operator++() {
74 while( sym.valid() && !(free()>0) ) { ++sym; }
79 void getFirst(OutEdgeIt& e, NodeIt v) const {
80 e=OutEdgeIt(*this, v);
82 void getFirst(EachNodeIt& v) const { G.getFirst(v); }
84 template< typename It >
91 template< typename It >
92 It first(NodeIt v) const {
98 NodeIt tail(EdgeIt e) const { return G.aNode(e.sym); }
99 NodeIt head(EdgeIt e) const { return G.bNode(e.sym); }
101 NodeIt aNode(OutEdgeIt e) const { return G.aNode(e.sym); }
102 NodeIt bNode(OutEdgeIt e) const { return G.bNode(e.sym); }
104 int id(NodeIt v) const { return G.id(v); }
106 template <typename S>
108 typename Graph::NodeMap<S> node_map;
110 NodeMap(const ResGraph<Graph, Number, FlowMap, CapacityMap>& _G) : node_map(_G.G) { }
111 NodeMap(const ResGraph<Graph, Number, FlowMap, CapacityMap>& _G, S a) : node_map(_G.G, a) { }
112 void set(NodeIt nit, S a) { node_map.set(nit, a); }
113 S get(NodeIt nit) const { return node_map.get(nit); }
114 S& operator[](NodeIt nit) { return node_map[nit]; }
115 const S& operator[](NodeIt nit) const { return node_map[nit]; }
121 template<typename Graph, typename Number, typename FlowMap, typename CapacityMap>
124 typedef typename Graph::NodeIt NodeIt;
125 typedef typename Graph::EachNodeIt EachNodeIt;
127 //typedef typename Graph::SymEdgeIt OldSymEdgeIt;
128 typedef typename Graph::OutEdgeIt OldOutEdgeIt;
129 typedef typename Graph::InEdgeIt OldInEdgeIt;
133 const CapacityMap& capacity;
135 ResGraph2(const Graph& _G, FlowMap& _flow,
136 const CapacityMap& _capacity) :
137 G(_G), flow(_flow), capacity(_capacity) { }
142 friend class OutEdgeIt;
145 friend class ResGraph2<Graph, Number, FlowMap, CapacityMap>;
147 const ResGraph2<Graph, Number, FlowMap, CapacityMap>* resG;
151 bool out_or_in; //true, iff out
153 EdgeIt() : out_or_in(true) { }
154 Number free() const {
156 return (resG->capacity.get(out)-resG->flow.get(out));
158 return (resG->flow.get(in));
162 return out_or_in && out.valid() || in.valid(); }
163 void augment(Number a) const {
165 resG->flow.set(out, resG->flow.get(out)+a);
167 resG->flow.set(in, resG->flow.get(in)-a);
172 class OutEdgeIt : public EdgeIt {
173 friend class ResGraph2<Graph, Number, FlowMap, CapacityMap>;
177 OutEdgeIt(const ResGraph2<Graph, Number, FlowMap, CapacityMap>& _resG, NodeIt v) {
179 out=resG->G.template first<OldOutEdgeIt>(v);
180 while( out.valid() && !(free()>0) ) { ++out; }
183 in=resG->G.template first<OldInEdgeIt>(v);
184 while( in.valid() && !(free()>0) ) { ++in; }
188 OutEdgeIt& operator++() {
190 NodeIt v=resG->G.aNode(out);
192 while( out.valid() && !(free()>0) ) { ++out; }
195 in=resG->G.template first<OldInEdgeIt>(v);
196 while( in.valid() && !(free()>0) ) { ++in; }
200 while( in.valid() && !(free()>0) ) { ++in; }
206 void getFirst(OutEdgeIt& e, NodeIt v) const {
207 e=OutEdgeIt(*this, v);
209 void getFirst(EachNodeIt& v) const { G.getFirst(v); }
211 template< typename It >
218 template< typename It >
219 It first(NodeIt v) const {
225 NodeIt tail(EdgeIt e) const {
226 return ((e.out_or_in) ? G.aNode(e.out) : G.aNode(e.in)); }
227 NodeIt head(EdgeIt e) const {
228 return ((e.out_or_in) ? G.bNode(e.out) : G.bNode(e.in)); }
230 NodeIt aNode(OutEdgeIt e) const {
231 return ((e.out_or_in) ? G.aNode(e.out) : G.aNode(e.in)); }
232 NodeIt bNode(OutEdgeIt e) const {
233 return ((e.out_or_in) ? G.bNode(e.out) : G.bNode(e.in)); }
235 int id(NodeIt v) const { return G.id(v); }
237 template <typename S>
239 typename Graph::NodeMap<S> node_map;
241 NodeMap(const ResGraph2<Graph, Number, FlowMap, CapacityMap>& _G) : node_map(_G.G) { }
242 NodeMap(const ResGraph2<Graph, Number, FlowMap, CapacityMap>& _G, S a) : node_map(_G.G, a) { }
243 void set(NodeIt nit, S a) { node_map.set(nit, a); }
244 S get(NodeIt nit) const { return node_map.get(nit); }
249 template <typename Graph, typename Number, typename FlowMap, typename CapacityMap>
252 typedef typename Graph::NodeIt NodeIt;
253 typedef typename Graph::EdgeIt EdgeIt;
254 typedef typename Graph::EachEdgeIt EachEdgeIt;
255 typedef typename Graph::OutEdgeIt OutEdgeIt;
256 typedef typename Graph::InEdgeIt InEdgeIt;
263 const CapacityMap* capacity;
264 typedef ResGraphWrapper<Graph, Number, FlowMap, CapacityMap > AugGraph;
265 typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
266 typedef typename AugGraph::EdgeIt AugEdgeIt;
268 //AugGraph res_graph;
269 //typedef typename AugGraph::NodeMap<bool> ReachedMap;
270 //typename AugGraph::NodeMap<AugEdgeIt> pred;
271 //typename AugGraph::NodeMap<Number> free;
273 MaxFlow(const Graph& _G, NodeIt _s, NodeIt _t, FlowMap& _flow, const CapacityMap& _capacity) :
274 G(&_G), s(_s), t(_t), flow(&_flow), capacity(&_capacity) //,
275 //res_graph(G, flow, capacity), pred(res_graph), free(res_graph)
277 bool augmentOnShortestPath() {
278 AugGraph res_graph(*G, *flow, *capacity);
281 typedef typename AugGraph::NodeMap<bool> ReachedMap;
282 BfsIterator5< AugGraph, /*AugOutEdgeIt,*/ ReachedMap > res_bfs(res_graph);
283 res_bfs.pushAndSetReached(s);
285 typename AugGraph::NodeMap<AugEdgeIt> pred(res_graph);
286 //filled up with invalid iterators
287 //pred.set(s, AugEdgeIt());
289 typename AugGraph::NodeMap<Number> free(res_graph);
291 //searching for augmenting path
292 while ( !res_bfs.finished() ) {
293 AugOutEdgeIt e=/*AugOutEdgeIt*/(res_bfs);
294 if (res_graph.valid(e) && res_bfs.isBNodeNewlyReached()) {
295 NodeIt v=res_graph.tail(e);
296 NodeIt w=res_graph.head(e);
298 if (res_graph.valid(pred.get(v))) {
299 free.set(w, std::min(free.get(v), res_graph.free(e)));
301 free.set(w, res_graph.free(e));
303 if (res_graph.head(e)==t) { _augment=true; break; }
307 } //end of searching augmenting path
311 Number augment_value=free.get(t);
312 while (res_graph.valid(pred.get(n))) {
313 AugEdgeIt e=pred.get(n);
314 res_graph.augment(e, augment_value);
315 //e.augment(augment_value);
323 template<typename MutableGraph> bool augmentOnBlockingFlow() {
326 AugGraph res_graph(*G, *flow, *capacity);
328 typedef typename AugGraph::NodeMap<bool> ReachedMap;
329 BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > bfs(res_graph);
331 bfs.pushAndSetReached(s);
332 typename AugGraph::NodeMap<int> dist(res_graph); //filled up with 0's
333 while ( !bfs.finished() ) {
334 AugOutEdgeIt e=/*AugOutEdgeIt*/(bfs);
335 if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
336 dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
340 } //computing distances from s in the residual graph
343 typename AugGraph::NodeMap<NodeIt> res_graph_to_F(res_graph);
344 for(typename AugGraph::EachNodeIt n=res_graph.template first<typename AugGraph::EachNodeIt>(); res_graph.valid(n); res_graph.next(n)) {
345 res_graph_to_F.set(n, F.addNode());
348 typename MutableGraph::NodeIt sF=res_graph_to_F.get(s);
349 typename MutableGraph::NodeIt tF=res_graph_to_F.get(t);
351 typename MutableGraph::EdgeMap<AugEdgeIt> original_edge(F);
352 typename MutableGraph::EdgeMap<Number> residual_capacity(F);
354 //Making F to the graph containing the edges of the residual graph
355 //which are in some shortest paths
356 for(typename AugGraph::EachEdgeIt e=res_graph.template first<typename AugGraph::EachEdgeIt>(); res_graph.valid(e); res_graph.next(e)) {
357 if (dist.get(res_graph.head(e))==dist.get(res_graph.tail(e))+1) {
358 typename MutableGraph::EdgeIt f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
359 original_edge.update();
360 original_edge.set(f, e);
361 residual_capacity.update();
362 residual_capacity.set(f, res_graph.free(e));
370 //computing blocking flow with dfs
371 typedef typename MutableGraph::NodeMap<bool> BlockingReachedMap;
372 DfsIterator4< MutableGraph, typename MutableGraph::OutEdgeIt, BlockingReachedMap > dfs(F);
373 typename MutableGraph::NodeMap<EdgeIt> pred(F); //invalid iterators
374 typename MutableGraph::NodeMap<Number> free(F);
376 dfs.pushAndSetReached(sF);
377 while (!dfs.finished()) {
379 if (F.valid(typename MutableGraph::OutEdgeIt(dfs))) {
380 if (dfs.isBNodeNewlyReached()) {
381 // std::cout << "OutEdgeIt: " << dfs;
382 // std::cout << " aNode: " << F.aNode(dfs);
383 // std::cout << " bNode: " << F.bNode(dfs) << " ";
385 typename MutableGraph::NodeIt v=F.aNode(dfs);
386 typename MutableGraph::NodeIt w=F.bNode(dfs);
388 if (F.valid(pred.get(v))) {
389 free.set(w, std::min(free.get(v), residual_capacity.get(dfs)));
391 free.set(w, residual_capacity.get(dfs));
394 //std::cout << "AUGMENTATION"<<std::endl;
401 F.erase(typename MutableGraph::OutEdgeIt(dfs));
407 typename MutableGraph::NodeIt n=tF;
408 Number augment_value=free.get(tF);
409 while (F.valid(pred.get(n))) {
410 typename MutableGraph::EdgeIt e=pred.get(n);
411 res_graph.augment(original_edge.get(e), augment_value);
412 //original_edge.get(e).augment(augment_value);
414 if (residual_capacity.get(e)==augment_value)
417 residual_capacity.set(e, residual_capacity.get(e)-augment_value);
425 bool augmentOnBlockingFlow2() {
428 //typedef ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> EAugGraph;
429 typedef FilterGraphWrapper< ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> > EAugGraph;
430 typedef typename EAugGraph::OutEdgeIt EAugOutEdgeIt;
431 typedef typename EAugGraph::EdgeIt EAugEdgeIt;
433 EAugGraph res_graph(*G, *flow, *capacity);
435 //std::cout << "meg jo1" << std::endl;
437 //typedef typename EAugGraph::NodeMap<bool> ReachedMap;
439 ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>,
440 ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt,
441 ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::NodeMap<bool> > bfs(res_graph);
443 //std::cout << "meg jo2" << std::endl;
445 bfs.pushAndSetReached(s);
446 //std::cout << "meg jo2.5" << std::endl;
448 //typename EAugGraph::NodeMap<int> dist(res_graph); //filled up with 0's
449 typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::
450 NodeMap<int>& dist=res_graph.dist;
451 //std::cout << "meg jo2.6" << std::endl;
453 while ( !bfs.finished() ) {
454 ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt e=bfs;
455 // EAugOutEdgeIt e=/*AugOutEdgeIt*/(bfs);
456 //if (res_graph.valid(e)) {
457 // std::cout<<"a:"<<res_graph.tail(e)<<"b:"<<res_graph.head(e)<<std::endl;
459 if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
460 dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
464 } //computing distances from s in the residual graph
467 //std::cout << "meg jo3" << std::endl;
469 // typedef typename EAugGraph::EachNodeIt EAugEachNodeIt;
470 // for(EAugEachNodeIt n=res_graph.template first<EAugEachNodeIt>(); res_graph.valid(n); res_graph.next(n)) {
471 // std::cout << "dist: " << dist.get(n) << std::endl;
477 // std::cout << "new iteration"<< std::endl;
480 //computing blocking flow with dfs
481 typedef typename EAugGraph::NodeMap<bool> BlockingReachedMap;
482 DfsIterator4< EAugGraph, EAugOutEdgeIt, BlockingReachedMap >
484 typename EAugGraph::NodeMap<EAugEdgeIt> pred(res_graph); //invalid iterators
485 typename EAugGraph::NodeMap<Number> free(res_graph);
487 dfs.pushAndSetReached(s);
488 while (!dfs.finished()) {
490 if (res_graph.valid(EAugOutEdgeIt(dfs))) {
491 if (dfs.isBNodeNewlyReached()) {
492 // std::cout << "OutEdgeIt: " << dfs;
493 // std::cout << " aNode: " << res_graph.aNode(dfs);
494 // std::cout << " res cap: " << EAugOutEdgeIt(dfs).free();
495 // std::cout << " bNode: " << res_graph.bNode(dfs) << " ";
497 typename EAugGraph::NodeIt v=res_graph.aNode(dfs);
498 typename EAugGraph::NodeIt w=res_graph.bNode(dfs);
500 pred.set(w, EAugOutEdgeIt(dfs));
502 //std::cout << EAugOutEdgeIt(dfs).free() << std::endl;
503 if (res_graph.valid(pred.get(v))) {
504 free.set(w, std::min(free.get(v), res_graph.free(/*EAugOutEdgeIt*/(dfs))));
506 free.set(w, res_graph.free(/*EAugOutEdgeIt*/(dfs)));
510 // std::cout << "t is reached, AUGMENTATION"<<std::endl;
516 // std::cout << "<<DELETE ";
517 // std::cout << " aNode: " << res_graph.aNode(dfs);
518 // std::cout << " res cap: " << EAugOutEdgeIt(dfs).free();
519 // std::cout << " bNode: " << res_graph.bNode(dfs) << " ";
520 // std::cout << "DELETE>> ";
522 res_graph.erase(dfs);
529 typename EAugGraph::NodeIt n=t;
530 Number augment_value=free.get(t);
531 // std::cout << "av:" << augment_value << std::endl;
532 while (res_graph.valid(pred.get(n))) {
533 EAugEdgeIt e=pred.get(n);
534 res_graph.augment(e, augment_value);
535 //e.augment(augment_value);
537 if (res_graph.free(e)==0)
547 //int num_of_augmentations=0;
548 while (augmentOnShortestPath()) {
549 //while (augmentOnBlockingFlow<MutableGraph>()) {
550 //std::cout << ++num_of_augmentations << " ";
551 //std::cout<<std::endl;
554 template<typename MutableGraph> void run() {
555 //int num_of_augmentations=0;
556 //while (augmentOnShortestPath()) {
557 while (augmentOnBlockingFlow<MutableGraph>()) {
558 //std::cout << ++num_of_augmentations << " ";
559 //std::cout<<std::endl;
565 for(G->getFirst(e, s); G->valid(e); G->next(e)) {
573 // template <typename Graph, typename Number, typename FlowMap, typename CapacityMap>
576 // typedef typename Graph::NodeIt NodeIt;
577 // typedef typename Graph::EdgeIt EdgeIt;
578 // typedef typename Graph::EachEdgeIt EachEdgeIt;
579 // typedef typename Graph::OutEdgeIt OutEdgeIt;
580 // typedef typename Graph::InEdgeIt InEdgeIt;
583 // std::list<NodeIt>& S;
584 // std::list<NodeIt>& T;
586 // const CapacityMap& capacity;
587 // typedef ResGraphWrapper<Graph, Number, FlowMap, CapacityMap > AugGraph;
588 // typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
589 // typedef typename AugGraph::EdgeIt AugEdgeIt;
590 // typename Graph::NodeMap<bool> SMap;
591 // typename Graph::NodeMap<bool> TMap;
593 // MaxFlow2(const Graph& _G, std::list<NodeIt>& _S, std::list<NodeIt>& _T, FlowMap& _flow, const CapacityMap& _capacity) : G(_G), S(_S), T(_T), flow(_flow), capacity(_capacity), SMap(_G), TMap(_G) {
594 // for(typename std::list<NodeIt>::const_iterator i=S.begin();
595 // i!=S.end(); ++i) {
596 // SMap.set(*i, true);
598 // for (typename std::list<NodeIt>::const_iterator i=T.begin();
599 // i!=T.end(); ++i) {
600 // TMap.set(*i, true);
604 // AugGraph res_graph(G, flow, capacity);
605 // bool _augment=false;
606 // NodeIt reached_t_node;
608 // typedef typename AugGraph::NodeMap<bool> ReachedMap;
609 // BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > res_bfs(res_graph);
610 // for(typename std::list<NodeIt>::const_iterator i=S.begin();
611 // i!=S.end(); ++i) {
612 // res_bfs.pushAndSetReached(*i);
614 // //res_bfs.pushAndSetReached(s);
616 // typename AugGraph::NodeMap<AugEdgeIt> pred(res_graph);
617 // //filled up with invalid iterators
619 // typename AugGraph::NodeMap<Number> free(res_graph);
621 // //searching for augmenting path
622 // while ( !res_bfs.finished() ) {
623 // AugOutEdgeIt e=/*AugOutEdgeIt*/(res_bfs);
624 // if (e.valid() && res_bfs.isBNodeNewlyReached()) {
625 // NodeIt v=res_graph.tail(e);
626 // NodeIt w=res_graph.head(e);
628 // if (pred.get(v).valid()) {
629 // free.set(w, std::min(free.get(v), e.free()));
631 // free.set(w, e.free());
633 // if (TMap.get(res_graph.head(e))) {
635 // reached_t_node=res_graph.head(e);
641 // } //end of searching augmenting path
644 // NodeIt n=reached_t_node;
645 // Number augment_value=free.get(reached_t_node);
646 // while (pred.get(n).valid()) {
647 // AugEdgeIt e=pred.get(n);
648 // e.augment(augment_value);
649 // n=res_graph.tail(e);
656 // while (augment()) { }
658 // Number flowValue() {
660 // for(typename std::list<NodeIt>::const_iterator i=S.begin();
661 // i!=S.end(); ++i) {
662 // for(OutEdgeIt e=G.template first<OutEdgeIt>(*i); e.valid(); ++e) {
665 // for(InEdgeIt e=G.template first<InEdgeIt>(*i); e.valid(); ++e) {
677 #endif //EDMONDS_KARP_HH