src/work/edmonds_karp.hh
changeset 134 e606071614f0
parent 127 dcace15b1874
child 135 1e5060d1fa1d
equal deleted inserted replaced
8:1bb8dd30d045 9:f0d33e345494
     4 #include <algorithm>
     4 #include <algorithm>
     5 #include <list>
     5 #include <list>
     6 #include <iterator>
     6 #include <iterator>
     7 
     7 
     8 #include <bfs_iterator.hh>
     8 #include <bfs_iterator.hh>
     9 #include <time_measure.h>
     9 //#include <time_measure.h>
    10 
    10 
    11 namespace hugo {
    11 namespace hugo {
    12 
    12 
    13   template<typename Graph, typename Number, typename FlowMap, typename CapacityMap>
    13   template<typename Graph, typename Number, typename FlowMap, typename CapacityMap>
    14   class ResGraph {
    14   class ResGraph {
   279       OldOutEdgeIt out;
   279       OldOutEdgeIt out;
   280       OldInEdgeIt in;
   280       OldInEdgeIt in;
   281       bool out_or_in; //1, iff out
   281       bool out_or_in; //1, iff out
   282     public:
   282     public:
   283       EdgeIt() : out_or_in(1) { } 
   283       EdgeIt() : out_or_in(1) { } 
       
   284       //EdgeIt(const EdgeIt& e) : G(e.G), flow(e.flow), capacity(e.capacity), out(e.out), in(e.in), out_or_in(e.out_or_in) { }
   284       Number free() const { 
   285       Number free() const { 
   285 	if (out_or_in) { 
   286 	if (out_or_in) { 
   286 	  return (/*resG->*/capacity->get(out)-/*resG->*/flow->get(out)); 
   287 	  return (/*resG->*/capacity->get(out)-/*resG->*/flow->get(out)); 
   287 	} else { 
   288 	} else { 
   288 	  return (/*resG->*/flow->get(in)); 
   289 	  return (/*resG->*/flow->get(in)); 
   295 	  /*resG->*/flow->set(out, /*resG->*/flow->get(out)+a);
   296 	  /*resG->*/flow->set(out, /*resG->*/flow->get(out)+a);
   296 	} else { 
   297 	} else { 
   297 	  /*resG->*/flow->set(in, /*resG->*/flow->get(in)-a);
   298 	  /*resG->*/flow->set(in, /*resG->*/flow->get(in)-a);
   298 	}
   299 	}
   299       }
   300       }
       
   301       void print() { 
       
   302 	if (out_or_in) {
       
   303 	  std::cout << "out "; 
       
   304 	  if (out.valid()) 
       
   305 	    std::cout << G->id(G->tail(out)) << "--"<< G->id(out) <<"->"<< G->id(G->head(out)); 
       
   306 	  else 
       
   307 	    std::cout << "invalid"; 
       
   308 	}
       
   309 	else {
       
   310 	  std::cout << "in "; 
       
   311 	  if (in.valid()) 
       
   312 	    std::cout << G->id(G->head(in)) << "<-"<< G->id(in) <<"--"<< G->id(G->tail(in)); 
       
   313 	  else 
       
   314 	    std::cout << "invalid"; 
       
   315 	}
       
   316 	std::cout << std::endl;
       
   317       }
   300     };
   318     };
   301 
   319 
   302     class OutEdgeIt : public EdgeIt {
   320     class OutEdgeIt : public EdgeIt {
   303       friend class ResGraph3<Graph, Number, FlowMap, CapacityMap>;
   321       friend class ResGraph3<Graph, Number, FlowMap, CapacityMap>;
   304     public:
   322     public:
   306     private:
   324     private:
   307       OutEdgeIt(const Graph& _G, NodeIt v, FlowMap& _flow, const CapacityMap& _capacity) { 
   325       OutEdgeIt(const Graph& _G, NodeIt v, FlowMap& _flow, const CapacityMap& _capacity) { 
   308       	G=&_G;
   326       	G=&_G;
   309 	flow=&_flow;
   327 	flow=&_flow;
   310 	capacity=&_capacity;
   328 	capacity=&_capacity;
   311 	out=/*resG->*/G->template first<OldOutEdgeIt>(v);
   329 	//out=/*resG->*/G->template first<OldOutEdgeIt>(v);
       
   330 	G->getFirst(out, v);
   312 	while( out.valid() && !(free()>0) ) { ++out; }
   331 	while( out.valid() && !(free()>0) ) { ++out; }
   313 	if (!out.valid()) {
   332 	if (!out.valid()) {
   314 	  out_or_in=0;
   333 	  out_or_in=0;
   315 	  in=/*resG->*/G->template first<OldInEdgeIt>(v);
   334 	  //in=/*resG->*/G->template first<OldInEdgeIt>(v);
       
   335 	  G->getFirst(in, v);
   316 	  while( in.valid() && !(free()>0) ) { ++in; }
   336 	  while( in.valid() && !(free()>0) ) { ++in; }
   317 	}
   337 	}
   318       }
   338       }
   319     public:
   339     public:
   320       OutEdgeIt& operator++() { 
   340       OutEdgeIt& operator++() { 
   322 	  NodeIt v=/*resG->*/G->aNode(out);
   342 	  NodeIt v=/*resG->*/G->aNode(out);
   323 	  ++out;
   343 	  ++out;
   324 	  while( out.valid() && !(free()>0) ) { ++out; }
   344 	  while( out.valid() && !(free()>0) ) { ++out; }
   325 	  if (!out.valid()) {
   345 	  if (!out.valid()) {
   326 	    out_or_in=0;
   346 	    out_or_in=0;
   327 	    in=/*resG->*/G->template first<OldInEdgeIt>(v);
   347 	    G->getFirst(in, v); //=/*resG->*/G->template first<OldInEdgeIt>(v);
   328 	    while( in.valid() && !(free()>0) ) { ++in; }
   348 	    while( in.valid() && !(free()>0) ) { ++in; }
   329 	  }
   349 	  }
   330 	} else {
   350 	} else {
   331 	  ++in;
   351 	  ++in;
   332 	  while( in.valid() && !(free()>0) ) { ++in; } 
   352 	  while( in.valid() && !(free()>0) ) { ++in; } 
   333 	}
   353 	}
   334 	return *this; 
   354 	return *this; 
   335       }
   355       }
   336     };
   356     };
   337 
   357 
       
   358     class EachEdgeIt : public EdgeIt {
       
   359       typename Graph::EachNodeIt v;
       
   360     public:
       
   361       EachEdgeIt() { }
       
   362       //EachEdgeIt(const EachEdgeIt& e) : EdgeIt(e), v(e.v) { }
       
   363       EachEdgeIt(const Graph& _G, FlowMap& _flow, const CapacityMap& _capacity) { 
       
   364 	G=&_G;
       
   365 	flow=&_flow;
       
   366 	capacity=&_capacity;
       
   367 	out_or_in=1;
       
   368 	G->getFirst(v);
       
   369 	if (v.valid()) G->getFirst(out, v); else out=OldOutEdgeIt();
       
   370 	while (out.valid() && !(free()>0) ) { ++out; }
       
   371 	while (v.valid() && !out.valid()) { 
       
   372 	  ++v; 
       
   373 	  if (v.valid()) G->getFirst(out, v); 
       
   374 	  while (out.valid() && !(free()>0) ) { ++out; }
       
   375 	}
       
   376 	if (!out.valid()) {
       
   377 	  out_or_in=0;
       
   378 	  G->getFirst(v);
       
   379 	  if (v.valid()) G->getFirst(in, v); else in=OldInEdgeIt();
       
   380 	  while (in.valid() && !(free()>0) ) { ++in; }
       
   381 	  while (v.valid() && !in.valid()) { 
       
   382 	    ++v; 
       
   383 	    if (v.valid()) G->getFirst(in, v); 
       
   384 	    while (in.valid() && !(free()>0) ) { ++in; }
       
   385 	  }
       
   386 	}
       
   387       }
       
   388       EachEdgeIt& operator++() { 
       
   389 	if (out_or_in) {
       
   390 	  ++out;
       
   391 	  while (out.valid() && !(free()>0) ) { ++out; }
       
   392 	  while (v.valid() && !out.valid()) { 
       
   393 	    ++v; 
       
   394 	    if (v.valid()) G->getFirst(out, v); 
       
   395 	    while (out.valid() && !(free()>0) ) { ++out; }
       
   396 	  }
       
   397 	  if (!out.valid()) {
       
   398 	    out_or_in=0;
       
   399 	    G->getFirst(v);
       
   400 	    if (v.valid()) G->getFirst(in, v); else in=OldInEdgeIt();
       
   401 	    while (in.valid() && !(free()>0) ) { ++in; }
       
   402 	    while (v.valid() && !in.valid()) { 
       
   403 	      ++v; 
       
   404 	      if (v.valid()) G->getFirst(in, v); 
       
   405 	      while (in.valid() && !(free()>0) ) { ++in; }
       
   406 	    }  
       
   407 	  }
       
   408 	} else {
       
   409 	  ++in;
       
   410 	  while (in.valid() && !(free()>0) ) { ++in; }
       
   411 	  while (v.valid() && !in.valid()) { 
       
   412 	    ++v; 
       
   413 	    if (v.valid()) G->getFirst(in, v); 
       
   414 	    while (in.valid() && !(free()>0) ) { ++in; }
       
   415 	  }
       
   416 	}
       
   417 	return *this;
       
   418       }
       
   419     };
       
   420 
   338     void getFirst(OutEdgeIt& e, NodeIt v) const { 
   421     void getFirst(OutEdgeIt& e, NodeIt v) const { 
   339       e=OutEdgeIt(G, v, flow, capacity); 
   422       e=OutEdgeIt(G, v, flow, capacity); 
       
   423     }
       
   424     void getFirst(EachEdgeIt& e) const { 
       
   425       e=EachEdgeIt(G, flow, capacity); 
   340     }
   426     }
   341     void getFirst(EachNodeIt& v) const { G.getFirst(v); }
   427     void getFirst(EachNodeIt& v) const { G.getFirst(v); }
   342     
   428     
   343     template< typename It >
   429     template< typename It >
   344     It first() const { 
   430     It first() const { 
   399     typedef typename AugGraph::EdgeIt AugEdgeIt;
   485     typedef typename AugGraph::EdgeIt AugEdgeIt;
   400 
   486 
   401     //AugGraph res_graph;    
   487     //AugGraph res_graph;    
   402     //typedef typename AugGraph::NodeMap<bool> ReachedMap;
   488     //typedef typename AugGraph::NodeMap<bool> ReachedMap;
   403     //typename AugGraph::NodeMap<AugEdgeIt> pred; 
   489     //typename AugGraph::NodeMap<AugEdgeIt> pred; 
   404     //typename AugGraph::NodeMap<int> free;
   490     //typename AugGraph::NodeMap<Number> free;
   405   public:
   491   public:
   406     MaxFlow(const Graph& _G, NodeIt _s, NodeIt _t, FlowMap& _flow, const CapacityMap& _capacity) : 
   492     MaxFlow(const Graph& _G, NodeIt _s, NodeIt _t, FlowMap& _flow, const CapacityMap& _capacity) : 
   407       G(_G), s(_s), t(_t), flow(_flow), capacity(_capacity) //,  
   493       G(_G), s(_s), t(_t), flow(_flow), capacity(_capacity) //,  
   408       //res_graph(G, flow, capacity), pred(res_graph), free(res_graph) 
   494       //res_graph(G, flow, capacity), pred(res_graph), free(res_graph) 
   409       { }
   495       { }
   410     bool augment() {
   496     bool augmentOnShortestPath() {
   411       AugGraph res_graph(G, flow, capacity);
   497       AugGraph res_graph(G, flow, capacity);
   412       bool _augment=false;
   498       bool _augment=false;
   413       
   499       
   414       typedef typename AugGraph::NodeMap<bool> ReachedMap;
   500       typedef typename AugGraph::NodeMap<bool> ReachedMap;
   415       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > res_bfs(res_graph);
   501       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > res_bfs(res_graph);
   417 	
   503 	
   418       typename AugGraph::NodeMap<AugEdgeIt> pred(res_graph); 
   504       typename AugGraph::NodeMap<AugEdgeIt> pred(res_graph); 
   419       //filled up with invalid iterators
   505       //filled up with invalid iterators
   420       //pred.set(s, AugEdgeIt());
   506       //pred.set(s, AugEdgeIt());
   421       
   507       
   422       typename AugGraph::NodeMap<int> free(res_graph);
   508       typename AugGraph::NodeMap<Number> free(res_graph);
   423 	
   509 	
   424       //searching for augmenting path
   510       //searching for augmenting path
   425       while ( !res_bfs.finished() ) { 
   511       while ( !res_bfs.finished() ) { 
   426 	AugOutEdgeIt e=AugOutEdgeIt(res_bfs);
   512 	AugOutEdgeIt e=AugOutEdgeIt(res_bfs);
   427 	if (e.valid() && res_bfs.isBNodeNewlyReached()) {
   513 	if (e.valid() && res_bfs.isBNodeNewlyReached()) {
   449 	}
   535 	}
   450       }
   536       }
   451 
   537 
   452       return _augment;
   538       return _augment;
   453     }
   539     }
   454     bool augmentWithBlockingFlow() {
   540     template<typename MutableGraph> bool augmentOnBlockingFlow() {
   455       BfsIterator4< Graph, OutEdgeIt, typename Graph::NodeMap<bool> > bfs(G);
   541       bool _augment=false;
       
   542 
       
   543       AugGraph res_graph(G, flow, capacity);
       
   544 
       
   545       typedef typename AugGraph::NodeMap<bool> ReachedMap;
       
   546       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > bfs(res_graph);
       
   547 
   456       bfs.pushAndSetReached(s);
   548       bfs.pushAndSetReached(s);
   457       typename Graph::NodeMap<int> dist(G); //filled up with 0's
   549       typename AugGraph::NodeMap<int> dist(res_graph); //filled up with 0's
   458       while ( !bfs.finished() ) { 
   550       while ( !bfs.finished() ) { 
   459 	OutEdgeIt e=OutEdgeIt(bfs);
   551 	AugOutEdgeIt e=AugOutEdgeIt(bfs);
   460 	if (e.valid() && bfs.isBNodeNewlyReached()) {
   552 	if (e.valid() && bfs.isBNodeNewlyReached()) {
   461 	  dist.set(G.head(e), dist.get(G.tail(e))+1);
   553 	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
   462 	  //NodeIt v=res_graph.tail(e);
       
   463 	  //NodeIt w=res_graph.head(e);
       
   464 	  //pred.set(w, e);
       
   465 	  //if (pred.get(v).valid()) {
       
   466 	  //  free.set(w, std::min(free.get(v), e.free()));
       
   467 	  //} else {
       
   468 	  //  free.set(w, e.free()); 
       
   469 	  //}
       
   470 	  //if (res_graph.head(e)==t) { _augment=true; break; }
       
   471 	}
   554 	}
   472 	
   555 	
   473 	++bfs;
   556 	++bfs;
   474       } //end of searching augmenting path
   557       } //computing distances from s in the residual graph
   475 
   558 
   476       double pre_time_copy=currTime();
       
   477       typedef Graph MutableGraph;
       
   478       MutableGraph F;
   559       MutableGraph F;
   479       typename Graph::NodeMap<NodeIt> G_to_F(G);
   560       typename AugGraph::NodeMap<NodeIt> res_graph_to_F(res_graph);
   480       for(typename Graph::EachNodeIt n=G.template first<typename Graph::EachNodeIt>(); n.valid(); ++n) {
   561       for(typename AugGraph::EachNodeIt n=res_graph.template first<typename AugGraph::EachNodeIt>(); n.valid(); ++n) {
   481 	G_to_F.set(n, F.addNode());
   562 	res_graph_to_F.set(n, F.addNode());
   482       }
   563       }
   483       for(typename Graph::EachEdgeIt e=G.template first<typename Graph::EachEdgeIt>(); e.valid(); ++e) {
   564       
   484 	if (dist.get(G.head(e))==dist.get(G.tail(e))+1) {
   565       typename MutableGraph::NodeIt sF=res_graph_to_F.get(s);
   485 	  F.addEdge(G_to_F.get(G.tail(e)), G_to_F.get(G.head(e)));
   566       typename MutableGraph::NodeIt tF=res_graph_to_F.get(t);
   486 	}
   567 
   487       }
   568       typename MutableGraph::EdgeMap<AugEdgeIt> original_edge(F);
   488       double post_time_copy=currTime();
   569       typename MutableGraph::EdgeMap<Number> free_on_edge(F);
   489       std::cout << "copy time: " << post_time_copy-pre_time_copy << " sec"<< std::endl; 
   570 
   490 
   571       //Making F to the graph containing the edges of the residual graph 
   491       return 0;
   572       //which are in some shortest paths
       
   573       for(typename AugGraph::EachEdgeIt e=res_graph.template first<typename AugGraph::EachEdgeIt>(); e.valid(); ++e) {
       
   574 	if (dist.get(res_graph.head(e))==dist.get(res_graph.tail(e))+1) {
       
   575 	  typename MutableGraph::EdgeIt f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
       
   576 	  original_edge.update();
       
   577 	  original_edge.set(f, e);
       
   578 	  free_on_edge.update();
       
   579 	  free_on_edge.set(f, e.free());
       
   580 	} 
       
   581       }
       
   582 
       
   583       bool __augment=true;
       
   584 
       
   585       while (__augment) {
       
   586 	__augment=false;
       
   587 	//computing blocking flow with dfs
       
   588 	typedef typename MutableGraph::NodeMap<bool> BlockingReachedMap;
       
   589 	DfsIterator4< MutableGraph, typename MutableGraph::OutEdgeIt, BlockingReachedMap > dfs(F);
       
   590 	typename MutableGraph::NodeMap<EdgeIt> pred(F); //invalid iterators
       
   591 	typename MutableGraph::NodeMap<Number> free(F);
       
   592 
       
   593 	dfs.pushAndSetReached(sF);      
       
   594 	while (!dfs.finished()) {
       
   595 	  ++dfs;
       
   596 	  if (typename MutableGraph::OutEdgeIt(dfs).valid()) {
       
   597 	    //std::cout << "OutEdgeIt: " << dfs; 
       
   598 	    //std::cout << " aNode: " << F.aNode(dfs); 
       
   599 	    //std::cout << " bNode: " << F.bNode(dfs) << " ";
       
   600 	  
       
   601 	    typename MutableGraph::NodeIt v=F.aNode(dfs);
       
   602 	    typename MutableGraph::NodeIt w=F.bNode(dfs);
       
   603 	    pred.set(w, dfs);
       
   604 	    if (pred.get(v).valid()) {
       
   605 	      free.set(w, std::min(free.get(v), free_on_edge.get(dfs)));
       
   606 	    } else {
       
   607 	      free.set(w, free_on_edge.get(dfs)/*original_edge.get(dfs).free()*/); 
       
   608 	    }
       
   609 	    if (w==tF) { 
       
   610 	      //std::cout << "AUGMENTATION"<<std::endl;
       
   611 	      __augment=true; 
       
   612 	      _augment=true;
       
   613 	      break; 
       
   614 	    }
       
   615 	  } else { 
       
   616 	    //std::cout << "OutEdgeIt: " << "invalid"; 
       
   617 	    //std::cout << " aNode: " << dfs.aNode(); 
       
   618 	    //std::cout << " bNode: " << "invalid" << " ";
       
   619 	  }
       
   620 	  if (dfs.isBNodeNewlyReached()) { 
       
   621 	    //std::cout << "bNodeIsNewlyReached ";
       
   622 	  } else { 
       
   623 	    //std::cout << "bNodeIsNotNewlyReached ";
       
   624 	    if (typename MutableGraph::OutEdgeIt(dfs).valid()) {
       
   625 	      //std::cout << "DELETE ";
       
   626 	      F.erase(typename MutableGraph::OutEdgeIt(dfs));
       
   627 	    }
       
   628 	  } 
       
   629 	  //if (dfs.isANodeExamined()) { 
       
   630 	    //std::cout << "aNodeIsExamined ";
       
   631 	    //} else { 
       
   632 	    //std::cout << "aNodeIsNotExamined ";
       
   633 	    //} 
       
   634 	  //std::cout<<std::endl;
       
   635 	}
       
   636 
       
   637 	if (__augment) {
       
   638 	  typename MutableGraph::NodeIt n=tF;
       
   639 	  Number augment_value=free.get(tF);
       
   640 	  while (pred.get(n).valid()) { 
       
   641 	    typename MutableGraph::EdgeIt e=pred.get(n);
       
   642 	    original_edge.get(e).augment(augment_value); 
       
   643 	    n=F.tail(e);
       
   644 	    F.erase(e);
       
   645 	  }
       
   646 	}
       
   647 
       
   648       }
       
   649             
       
   650       return _augment;
   492     }
   651     }
   493     void run() {
   652     void run() {
   494       //int num_of_augmentations=0;
   653       //int num_of_augmentations=0;
   495       while (augment()) { 
   654       while (augmentOnShortestPath()) { 
       
   655 	//while (augmentOnBlockingFlow<MutableGraph>()) { 
       
   656 	//std::cout << ++num_of_augmentations << " ";
       
   657 	//std::cout<<std::endl;
       
   658       } 
       
   659     }
       
   660     template<typename MutableGraph> void run() {
       
   661       //int num_of_augmentations=0;
       
   662       //while (augmentOnShortestPath()) { 
       
   663 	while (augmentOnBlockingFlow<MutableGraph>()) { 
   496 	//std::cout << ++num_of_augmentations << " ";
   664 	//std::cout << ++num_of_augmentations << " ";
   497 	//std::cout<<std::endl;
   665 	//std::cout<<std::endl;
   498       } 
   666       } 
   499     }
   667     }
   500     Number flowValue() { 
   668     Number flowValue() { 
   597       //res_bfs.pushAndSetReached(s);
   765       //res_bfs.pushAndSetReached(s);
   598 	
   766 	
   599       typename AugGraph::NodeMap<AugEdgeIt> pred(res_graph); 
   767       typename AugGraph::NodeMap<AugEdgeIt> pred(res_graph); 
   600       //filled up with invalid iterators
   768       //filled up with invalid iterators
   601       
   769       
   602       typename AugGraph::NodeMap<int> free(res_graph);
   770       typename AugGraph::NodeMap<Number> free(res_graph);
   603 	
   771 	
   604       //searching for augmenting path
   772       //searching for augmenting path
   605       while ( !res_bfs.finished() ) { 
   773       while ( !res_bfs.finished() ) { 
   606 	AugOutEdgeIt e=AugOutEdgeIt(res_bfs);
   774 	AugOutEdgeIt e=AugOutEdgeIt(res_bfs);
   607 	if (e.valid() && res_bfs.isBNodeNewlyReached()) {
   775 	if (e.valid() && res_bfs.isBNodeNewlyReached()) {