lemon-0.x: src/work/marci/augmenting

     1 // -*- C++ -*-

     2 #ifndef HUGO_AUGMENTING_FLOW_H

     3 #define HUGO_AUGMENTING_FLOW_H

     5 #include <vector>

     6 #include <iostream>

     8 #include <hugo/graph_wrapper.h>

     9 #include <bfs_dfs.h>

    10 #include <hugo/invalid.h>

    11 #include <hugo/maps.h>

    12 #include <hugo/tight_edge_filter_map.h>

    14 /// \file

    15 /// \brief Maximum flow algorithms.

    16 /// \ingroup galgs

    18 namespace hugo {

    20   /// \addtogroup galgs

    21   /// @{

    22   /// Class for augmenting path flow algorithms.

    24   /// This class provides various algorithms for finding a flow of

    25   /// maximum value in a directed graph. The \e source node, the \e

    26   /// target node, the \e capacity of the edges and the \e starting \e

    27   /// flow value of the edges should be passed to the algorithm through the

    28   /// constructor.

    29 //   /// It is possible to change these quantities using the

    30 //   /// functions \ref resetSource, \ref resetTarget, \ref resetCap and

    31 //   /// \ref resetFlow. Before any subsequent runs of any algorithm of

    32 //   /// the class \ref resetFlow should be called.

    34   /// After running an algorithm of the class, the actual flow value

    35   /// can be obtained by calling \ref flowValue(). The minimum

    36   /// value cut can be written into a \c node map of \c bools by

    37   /// calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes

    38   /// the inclusionwise minimum and maximum of the minimum value

    39   /// cuts, resp.)

    40   ///\param Graph The directed graph type the algorithm runs on.

    41   ///\param Num The number type of the capacities and the flow values.

    42   ///\param CapMap The capacity map type.

    43   ///\param FlowMap The flow map type.

    44   ///\author Marton Makai

    45   template <typename Graph, typename Num,

    46 	    typename CapMap=typename Graph::template EdgeMap<Num>,

    47             typename FlowMap=typename Graph::template EdgeMap<Num> >

    48   class AugmentingFlow {

    49   protected:

    50     typedef typename Graph::Node Node;

    51     typedef typename Graph::NodeIt NodeIt;

    52     typedef typename Graph::EdgeIt EdgeIt;

    53     typedef typename Graph::OutEdgeIt OutEdgeIt;

    54     typedef typename Graph::InEdgeIt InEdgeIt;

    56     const Graph* g;

    57     Node s;

    58     Node t;

    59     const CapMap* capacity;

    60     FlowMap* flow;

    61 //    int n;      //the number of nodes of G

    62     typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;

    63     //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;

    64     typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;

    65     typedef typename ResGW::Edge ResGWEdge;

    66     //typedef typename ResGW::template NodeMap<bool> ReachedMap;

    67     typedef typename Graph::template NodeMap<int> ReachedMap;

    69     //level works as a bool map in augmenting path algorithms and is

    70     //used by bfs for storing reached information.  In preflow, it

    71     //shows the levels of nodes.

    72     ReachedMap level;

    74   public:

    75     ///Indicates the property of the starting flow.

    77     ///Indicates the property of the starting flow. The meanings are as follows:

    78     ///- \c ZERO_FLOW: constant zero flow

    79     ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to

    80     ///the sum of the out-flows in every node except the \e source and

    81     ///the \e target.

    82     ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at

    83     ///least the sum of the out-flows in every node except the \e source.

    84     ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be

    85     ///set to the constant zero flow in the beginning of the algorithm in this case.

    86     enum FlowEnum{

    87       ZERO_FLOW,

    88       GEN_FLOW,

    89       PRE_FLOW,

    90       NO_FLOW

    91     };

    93     enum StatusEnum {

    94       AFTER_NOTHING,

    95       AFTER_AUGMENTING,

    96       AFTER_FAST_AUGMENTING,

    97       AFTER_PRE_FLOW_PHASE_1,

    98       AFTER_PRE_FLOW_PHASE_2

    99     };

   101     /// Don not needle this flag only if necessary.

   102     StatusEnum status;

   103     int number_of_augmentations;

   106     template<typename IntMap>

   107     class TrickyReachedMap {

   108     protected:

   109       IntMap* map;

   110       int* number_of_augmentations;

   111     public:

   112       TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) :

   113 	map(&_map), number_of_augmentations(&_number_of_augmentations) { }

   114       void set(const Node& n, bool b) {

   115 	if (b)

   116 	  map->set(n, *number_of_augmentations);

   117 	else

   118 	  map->set(n, *number_of_augmentations-1);

   119       }

   120       bool operator[](const Node& n) const {

   121 	return (*map)[n]==*number_of_augmentations;

   122       }

   123     };

   125     AugmentingFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,

   126 		   FlowMap& _flow) :

   127       g(&_G), s(_s), t(_t), capacity(&_capacity),

   128       flow(&_flow), //n(_G.nodeNum()),

   129       level(_G), //excess(_G,0),

   130       status(AFTER_NOTHING), number_of_augmentations(0) { }

   132     /// Starting from a flow, this method searches for an augmenting path

   133     /// according to the Edmonds-Karp algorithm

   134     /// and augments the flow on if any.

   135     /// The return value shows if the augmentation was succesful.

   136     bool augmentOnShortestPath();

   137     bool augmentOnShortestPath2();

   139     /// Starting from a flow, this method searches for an augmenting blocking

   140     /// flow according to Dinits' algorithm and augments the flow on if any.

   141     /// The blocking flow is computed in a physically constructed

   142     /// residual graph of type \c Mutablegraph.

   143     /// The return value show sif the augmentation was succesful.

   144     template<typename MutableGraph> bool augmentOnBlockingFlow();

   146     /// The same as \c augmentOnBlockingFlow<MutableGraph> but the

   147     /// residual graph is not constructed physically.

   148     /// The return value shows if the augmentation was succesful.

   149     bool augmentOnBlockingFlow2();

   151     template<typename _CutMap>

   152     void actMinCut(_CutMap& M) const {

   153       NodeIt v;

   154       switch (status) {

   155 	case AFTER_PRE_FLOW_PHASE_1:

   156 //	std::cout << "AFTER_PRE_FLOW_PHASE_1" << std::endl;

   157 // 	for(g->first(v); g->valid(v); g->next(v)) {

   158 // 	  if (level[v] < n) {

   159 // 	    M.set(v, false);

   160 // 	  } else {

   161 // 	    M.set(v, true);

   162 // 	  }

   163 // 	}

   164 	break;

   165       case AFTER_PRE_FLOW_PHASE_2:

   166 //	std::cout << "AFTER_PRE_FLOW_PHASE_2" << std::endl;

   167 	break;

   168       case AFTER_NOTHING:

   169 //	std::cout << "AFTER_NOTHING" << std::endl;

   170 	minMinCut(M);

   171 	break;

   172       case AFTER_AUGMENTING:

   173 //	std::cout << "AFTER_AUGMENTING" << std::endl;

   174 	for(g->first(v); v!=INVALID; ++v) {

   175 	  if (level[v]) {

   176 	    M.set(v, true);

   177 	  } else {

   178 	    M.set(v, false);

   179 	  }

   180 	}

   181 	break;

   182       case AFTER_FAST_AUGMENTING:

   183 //	std::cout << "AFTER_FAST_AUGMENTING" << std::endl;

   184 	for(g->first(v); v!=INVALID; ++v) {

   185 	  if (level[v]==number_of_augmentations) {

   186 	    M.set(v, true);

   187 	  } else {

   188 	    M.set(v, false);

   189 	  }

   190 	}

   191 	break;

   192       }

   193     }

   195     template<typename _CutMap>

   196     void minMinCut(_CutMap& M) const {

   197       std::queue<Node> queue;

   199       M.set(s,true);

   200       queue.push(s);

   202       while (!queue.empty()) {

   203         Node w=queue.front();

   204 	queue.pop();

   206 	OutEdgeIt e;

   207 	for(g->first(e,w) ; e!=INVALID; ++e) {

   208 	  Node v=g->head(e);

   209 	  if (!M[v] && (*flow)[e] < (*capacity)[e] ) {

   210 	    queue.push(v);

   211 	    M.set(v, true);

   212 	  }

   213 	}

   215 	InEdgeIt f;

   216 	for(g->first(f,w) ; f!=INVALID; ++f) {

   217 	  Node v=g->tail(f);

   218 	  if (!M[v] && (*flow)[f] > 0 ) {

   219 	    queue.push(v);

   220 	    M.set(v, true);

   221 	  }

   222 	}

   223       }

   224     }

   226     template<typename _CutMap>

   227     void minMinCut2(_CutMap& M) const {

   228       ResGW res_graph(*g, *capacity, *flow);

   229       BfsIterator<ResGW, _CutMap> bfs(res_graph, M);

   230       bfs.pushAndSetReached(s);

   231       while (!bfs.finished()) ++bfs;

   232     }

   234     Num flowValue() const {

   235       Num a=0;

   236       for (InEdgeIt e(*g, t); e!=INVALID; ++e) a+=(*flow)[e];

   237       for (OutEdgeIt e(*g, t); e!=INVALID; ++e) a-=(*flow)[e];

   238       return a;

   239       //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan

   240     }

   242   };

   246   template <typename Graph, typename Num, typename CapMap, typename FlowMap>

   247   bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()

   248   {

   249     ResGW res_graph(*g, *capacity, *flow);

   250     typename ResGW::ResCap res_cap(res_graph);

   252     bool _augment=false;

   254     //ReachedMap level(res_graph);

   255     for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0);

   256     BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);

   257     bfs.pushAndSetReached(s);

   259     typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);

   260     pred.set(s, INVALID);

   262     typename ResGW::template NodeMap<Num> free(res_graph);

   264     //searching for augmenting path

   265     while ( !bfs.finished() ) {

   266       ResGWEdge e=bfs;

   267       if (e!=INVALID && bfs.isBNodeNewlyReached()) {

   268 	Node v=res_graph.tail(e);

   269 	Node w=res_graph.head(e);

   270 	pred.set(w, e);

   271 	if (pred[v]!=INVALID) {

   272 	  free.set(w, std::min(free[v], res_cap[e]));

   273 	} else {

   274 	  free.set(w, res_cap[e]);

   275 	}

   276 	if (res_graph.head(e)==t) { _augment=true; break; }

   277       }

   279       ++bfs;

   280     } //end of searching augmenting path

   282     if (_augment) {

   283       Node n=t;

   284       Num augment_value=free[t];

   285       while (pred[n]!=INVALID) {

   286 	ResGWEdge e=pred[n];

   287 	res_graph.augment(e, augment_value);

   288 	n=res_graph.tail(e);

   289       }

   290     }

   292     status=AFTER_AUGMENTING;

   293     return _augment;

   294   }

   296   template <typename Graph, typename Num, typename CapMap, typename FlowMap>

   297   bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath2()

   298   {

   299     ResGW res_graph(*g, *capacity, *flow);

   300     typename ResGW::ResCap res_cap(res_graph);

   302     bool _augment=false;

   304     if (status!=AFTER_FAST_AUGMENTING) {

   305       for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0);

   306       number_of_augmentations=1;

   307     } else {

   308       ++number_of_augmentations;

   309     }

   310     TrickyReachedMap<ReachedMap>

   311       tricky_reached_map(level, number_of_augmentations);

   312     //ReachedMap level(res_graph);

   313 //    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);

   314     BfsIterator<ResGW, TrickyReachedMap<ReachedMap> >

   315       bfs(res_graph, tricky_reached_map);

   316     bfs.pushAndSetReached(s);

   318     typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);

   319     pred.set(s, INVALID);

   321     typename ResGW::template NodeMap<Num> free(res_graph);

   323     //searching for augmenting path

   324     while ( !bfs.finished() ) {

   325       ResGWEdge e=bfs;

   326       if (e!=INVALID && bfs.isBNodeNewlyReached()) {

   327 	Node v=res_graph.tail(e);

   328 	Node w=res_graph.head(e);

   329 	pred.set(w, e);

   330 	if (pred[v]!=INVALID) {

   331 	  free.set(w, std::min(free[v], res_cap[e]));

   332 	} else {

   333 	  free.set(w, res_cap[e]);

   334 	}

   335 	if (res_graph.head(e)==t) { _augment=true; break; }

   336       }

   338       ++bfs;

   339     } //end of searching augmenting path

   341     if (_augment) {

   342       Node n=t;

   343       Num augment_value=free[t];

   344       while (pred[n]!=INVALID) {

   345 	ResGWEdge e=pred[n];

   346 	res_graph.augment(e, augment_value);

   347 	n=res_graph.tail(e);

   348       }

   349     }

   351     status=AFTER_FAST_AUGMENTING;

   352     return _augment;

   353   }

   356   template <typename Graph, typename Num, typename CapMap, typename FlowMap>

   357   template<typename MutableGraph>

   358   bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()

   359   {

   360     typedef MutableGraph MG;

   361     bool _augment=false;

   363     ResGW res_graph(*g, *capacity, *flow);

   364     typename ResGW::ResCap res_cap(res_graph);

   366     //bfs for distances on the residual graph

   367     //ReachedMap level(res_graph);

   368     for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0);

   369     BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);

   370     bfs.pushAndSetReached(s);

   371     typename ResGW::template NodeMap<int>

   372       dist(res_graph); //filled up with 0's

   374     //F will contain the physical copy of the residual graph

   375     //with the set of edges which are on shortest paths

   376     MG F;

   377     typename ResGW::template NodeMap<typename MG::Node>

   378       res_graph_to_F(res_graph);

   379     {

   380       typename ResGW::NodeIt n;

   381       for(res_graph.first(n); n!=INVALID; ++n)

   382 	res_graph_to_F.set(n, F.addNode());

   383     }

   385     typename MG::Node sF=res_graph_to_F[s];

   386     typename MG::Node tF=res_graph_to_F[t];

   387     typename MG::template EdgeMap<ResGWEdge> original_edge(F);

   388     typename MG::template EdgeMap<Num> residual_capacity(F);

   390     while ( !bfs.finished() ) {

   391       ResGWEdge e=bfs;

   392       if (e!=INVALID) {

   393 	if (bfs.isBNodeNewlyReached()) {

   394 	  dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);

   395 	  typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],

   396 					res_graph_to_F[res_graph.head(e)]);

   397 	  //original_edge.update();

   398 	  original_edge.set(f, e);

   399 	  //residual_capacity.update();

   400 	  residual_capacity.set(f, res_cap[e]);

   401 	} else {

   402 	  if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {

   403 	    typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],

   404 					  res_graph_to_F[res_graph.head(e)]);

   405 	    //original_edge.update();

   406 	    original_edge.set(f, e);

   407 	    //residual_capacity.update();

   408 	    residual_capacity.set(f, res_cap[e]);

   409 	  }

   410 	}

   411       }

   412       ++bfs;

   413     } //computing distances from s in the residual graph

   415     bool __augment=true;

   417     while (__augment) {

   418       __augment=false;

   419       //computing blocking flow with dfs

   420       DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);

   421       typename MG::template NodeMap<typename MG::Edge> pred(F);

   422       pred.set(sF, INVALID);

   423       //invalid iterators for sources

   425       typename MG::template NodeMap<Num> free(F);

   427       dfs.pushAndSetReached(sF);

   428       while (!dfs.finished()) {

   429 	++dfs;

   430 	if (typename MG::Edge(dfs)!=INVALID) {

   431 	  if (dfs.isBNodeNewlyReached()) {

   432 	    typename MG::Node v=F.tail(dfs);

   433 	    typename MG::Node w=F.head(dfs);

   434 	    pred.set(w, dfs);

   435 	    if (pred[v]!=INVALID) {

   436 	      free.set(w, std::min(free[v], residual_capacity[dfs]));

   437 	    } else {

   438 	      free.set(w, residual_capacity[dfs]);

   439 	    }

   440 	    if (w==tF) {

   441 	      __augment=true;

   442 	      _augment=true;

   443 	      break;

   444 	    }

   446 	  } else {

   447 	    F.erase(typename MG::Edge(dfs));

   448 	  }

   449 	}

   450       }

   452       if (__augment) {

   453 	typename MG::Node n=tF;

   454 	Num augment_value=free[tF];

   455 	while (pred[n]!=INVALID) {

   456 	  typename MG::Edge e=pred[n];

   457 	  res_graph.augment(original_edge[e], augment_value);

   458 	  n=F.tail(e);

   459 	  if (residual_capacity[e]==augment_value)

   460 	    F.erase(e);

   461 	  else

   462 	    residual_capacity.set(e, residual_capacity[e]-augment_value);

   463 	}

   464       }

   466     }

   468     status=AFTER_AUGMENTING;

   469     return _augment;

   470   }

   472   /// Blocking flow augmentation without constructing the layered

   473   /// graph physically in which the blocking flow is computed.

   474   template <typename Graph, typename Num, typename CapMap, typename FlowMap>

   475   bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()

   476   {

   477     bool _augment=false;

   479     ResGW res_graph(*g, *capacity, *flow);

   480     typename ResGW::ResCap res_cap(res_graph);

   482     //Potential map, for distances from s

   483     typename ResGW::template NodeMap<int> potential(res_graph, 0);

   484     typedef ConstMap<typename ResGW::Edge, int> Const1Map;

   485     Const1Map const_1_map(1);

   486     TightEdgeFilterMap<ResGW, typename ResGW::template NodeMap<int>,

   487       Const1Map> tight_edge_filter(res_graph, potential, const_1_map);

   489     for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0);

   490     BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);

   491     bfs.pushAndSetReached(s);

   493     //computing distances from s in the residual graph

   494     while ( !bfs.finished() ) {

   495       ResGWEdge e=bfs;

   496       if (e!=INVALID && bfs.isBNodeNewlyReached())

   497 	potential.set(res_graph.head(e), potential[res_graph.tail(e)]+1);

   498       ++bfs;

   499     }

   501     //Subgraph containing the edges on some shortest paths

   502     //(i.e. tight edges)

   503     ConstMap<typename ResGW::Node, bool> true_map(true);

   504     typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,

   505       TightEdgeFilterMap<ResGW, typename ResGW::template NodeMap<int>,

   506       Const1Map> > FilterResGW;

   507     FilterResGW filter_res_graph(res_graph, true_map, tight_edge_filter);

   509     //Subgraph, which is able to delete edges which are already

   510     //met by the dfs

   511     typename FilterResGW::template NodeMap<typename FilterResGW::Edge>

   512       first_out_edges(filter_res_graph);

   513     for (typename FilterResGW::NodeIt v(filter_res_graph); v!=INVALID; ++v)

   514       first_out_edges.set

   515 	(v, typename FilterResGW::OutEdgeIt(filter_res_graph, v));

   517     typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::

   518       template NodeMap<typename FilterResGW::Edge> > ErasingResGW;

   519     ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);

   521     bool __augment=true;

   523     while (__augment) {

   525       __augment=false;

   526       //computing blocking flow with dfs

   527       DfsIterator< ErasingResGW,

   528 	typename ErasingResGW::template NodeMap<bool> >

   529 	dfs(erasing_res_graph);

   530       typename ErasingResGW::

   531 	template NodeMap<typename ErasingResGW::Edge> pred(erasing_res_graph);

   532       pred.set(s, INVALID);

   533       //invalid iterators for sources

   535       typename ErasingResGW::template NodeMap<Num>

   536 	free1(erasing_res_graph);

   538       dfs.pushAndSetReached

   539 	/// \bug hugo 0.2

   540 	(typename ErasingResGW::Node

   541 	 (typename FilterResGW::Node

   542 	  (typename ResGW::Node(s)

   543 	   )

   544 	  )

   545 	 );

   547       while (!dfs.finished()) {

   548 	++dfs;

   549 	if (typename ErasingResGW::Edge(dfs)!=INVALID) {

   550 	  if (dfs.isBNodeNewlyReached()) {

   552 	    typename ErasingResGW::Node v=erasing_res_graph.tail(dfs);

   553 	    typename ErasingResGW::Node w=erasing_res_graph.head(dfs);

   555 	    pred.set(w, typename ErasingResGW::Edge(dfs));

   556 	    if (pred[v]!=INVALID) {

   557 	      free1.set

   558 		(w, std::min(free1[v], res_cap

   559 			     [typename ErasingResGW::Edge(dfs)]));

   560 	    } else {

   561 	      free1.set

   562 		(w, res_cap

   563 		 [typename ErasingResGW::Edge(dfs)]);

   564 	    }

   566 	    if (w==t) {

   567 	      __augment=true;

   568 	      _augment=true;

   569 	      break;

   570 	    }

   571 	  } else {

   572 	    erasing_res_graph.erase(dfs);

   573 	  }

   574 	}

   575       }

   577       if (__augment) {

   578 	typename ErasingResGW::Node

   579 	  n=typename FilterResGW::Node(typename ResGW::Node(t));

   580 	Num augment_value=free1[n];

   581 	while (pred[n]!=INVALID) {

   582 	  typename ErasingResGW::Edge e=pred[n];

   583 	  res_graph.augment(e, augment_value);

   584 	  n=erasing_res_graph.tail(e);

   585 	  if (res_cap[e]==0)

   586 	    erasing_res_graph.erase(e);

   587 	}

   588       }

   590     } //while (__augment)

   592     status=AFTER_AUGMENTING;

   593     return _augment;

   594   }

   597 } //namespace hugo

   599 #endif //HUGO_AUGMENTING_FLOW_H

author	marci
	Wed, 22 Sep 2004 10:47:59 +0000
changeset 901	69a8e672acb1
parent 863	d27bbe17b0b8
child 921	818510fa3d99
permissions	-rw-r--r--