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 <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     bool _augment=false;

   252     //ReachedMap level(res_graph);

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

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

   255     bfs.pushAndSetReached(s);

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

   258     pred.set(s, INVALID);

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

   262     //searching for augmenting path

   263     while ( !bfs.finished() ) {

   264       ResGWEdge e=bfs;

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

   266 	Node v=res_graph.tail(e);

   267 	Node w=res_graph.head(e);

   268 	pred.set(w, e);

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

   270 	  free.set(w, std::min(free[v], res_graph.resCap(e)));

   271 	} else {

   272 	  free.set(w, res_graph.resCap(e));

   273 	}

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

   275       }

   277       ++bfs;

   278     } //end of searching augmenting path

   280     if (_augment) {

   281       Node n=t;

   282       Num augment_value=free[t];

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

   284 	ResGWEdge e=pred[n];

   285 	res_graph.augment(e, augment_value);

   286 	n=res_graph.tail(e);

   287       }

   288     }

   290     status=AFTER_AUGMENTING;

   291     return _augment;

   292   }

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

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

   296   {

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

   298     bool _augment=false;

   300     if (status!=AFTER_FAST_AUGMENTING) {

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

   302       number_of_augmentations=1;

   303     } else {

   304       ++number_of_augmentations;

   305     }

   306     TrickyReachedMap<ReachedMap>

   307       tricky_reached_map(level, number_of_augmentations);

   308     //ReachedMap level(res_graph);

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

   310     BfsIterator<ResGW, TrickyReachedMap<ReachedMap> >

   311       bfs(res_graph, tricky_reached_map);

   312     bfs.pushAndSetReached(s);

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

   315     pred.set(s, INVALID);

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

   319     //searching for augmenting path

   320     while ( !bfs.finished() ) {

   321       ResGWEdge e=bfs;

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

   323 	Node v=res_graph.tail(e);

   324 	Node w=res_graph.head(e);

   325 	pred.set(w, e);

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

   327 	  free.set(w, std::min(free[v], res_graph.resCap(e)));

   328 	} else {

   329 	  free.set(w, res_graph.resCap(e));

   330 	}

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

   332       }

   334       ++bfs;

   335     } //end of searching augmenting path

   337     if (_augment) {

   338       Node n=t;

   339       Num augment_value=free[t];

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

   341 	ResGWEdge e=pred[n];

   342 	res_graph.augment(e, augment_value);

   343 	n=res_graph.tail(e);

   344       }

   345     }

   347     status=AFTER_FAST_AUGMENTING;

   348     return _augment;

   349   }

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

   353   template<typename MutableGraph>

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

   355   {

   356     typedef MutableGraph MG;

   357     bool _augment=false;

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

   361     //bfs for distances on the residual graph

   362     //ReachedMap level(res_graph);

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

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

   365     bfs.pushAndSetReached(s);

   366     typename ResGW::template NodeMap<int>

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

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

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

   371     MG F;

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

   373       res_graph_to_F(res_graph);

   374     {

   375       typename ResGW::NodeIt n;

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

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

   378     }

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

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

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

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

   385     while ( !bfs.finished() ) {

   386       ResGWEdge e=bfs;

   387       if (e!=INVALID) {

   388 	if (bfs.isBNodeNewlyReached()) {

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

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

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

   392 	  //original_edge.update();

   393 	  original_edge.set(f, e);

   394 	  //residual_capacity.update();

   395 	  residual_capacity.set(f, res_graph.resCap(e));

   396 	} else {

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

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

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

   400 	    //original_edge.update();

   401 	    original_edge.set(f, e);

   402 	    //residual_capacity.update();

   403 	    residual_capacity.set(f, res_graph.resCap(e));

   404 	  }

   405 	}

   406       }

   407       ++bfs;

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

   410     bool __augment=true;

   412     while (__augment) {

   413       __augment=false;

   414       //computing blocking flow with dfs

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

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

   417       pred.set(sF, INVALID);

   418       //invalid iterators for sources

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

   422       dfs.pushAndSetReached(sF);

   423       while (!dfs.finished()) {

   424 	++dfs;

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

   426 	  if (dfs.isBNodeNewlyReached()) {

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

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

   429 	    pred.set(w, dfs);

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

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

   432 	    } else {

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

   434 	    }

   435 	    if (w==tF) {

   436 	      __augment=true;

   437 	      _augment=true;

   438 	      break;

   439 	    }

   441 	  } else {

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

   443 	  }

   444 	}

   445       }

   447       if (__augment) {

   448 	typename MG::Node n=tF;

   449 	Num augment_value=free[tF];

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

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

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

   453 	  n=F.tail(e);

   454 	  if (residual_capacity[e]==augment_value)

   455 	    F.erase(e);

   456 	  else

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

   458 	}

   459       }

   461     }

   463     status=AFTER_AUGMENTING;

   464     return _augment;

   465   }

   467   /// Blocking flow augmentation without constructing the layered

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

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

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

   471   {

   472     bool _augment=false;

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

   476     //Potential map, for distances from s

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

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

   479     Const1Map const_1_map(1);

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

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

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

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

   485     bfs.pushAndSetReached(s);

   487     //computing distances from s in the residual graph

   488     while ( !bfs.finished() ) {

   489       ResGWEdge e=bfs;

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

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

   492       ++bfs;

   493     }

   495     //Subgraph containing the edges on some shortest paths

   496     //(i.e. tight edges)

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

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

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

   500       Const1Map> > FilterResGW;

   501     FilterResGW filter_res_graph(res_graph, true_map, tight_edge_filter);

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

   504     //met by the dfs

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

   506       first_out_edges(filter_res_graph);

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

   508       first_out_edges.set

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

   511     typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::

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

   513     ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);

   515     bool __augment=true;

   517     while (__augment) {

   519       __augment=false;

   520       //computing blocking flow with dfs

   521       DfsIterator< ErasingResGW,

   522 	typename ErasingResGW::template NodeMap<bool> >

   523 	dfs(erasing_res_graph);

   524       typename ErasingResGW::

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

   526       pred.set(s, INVALID);

   527       //invalid iterators for sources

   529       typename ErasingResGW::template NodeMap<Num>

   530 	free1(erasing_res_graph);

   532       dfs.pushAndSetReached

   533 	/// \bug hugo 0.2

   534 	(typename ErasingResGW::Node

   535 	 (typename FilterResGW::Node

   536 	  (typename ResGW::Node(s)

   537 	   )

   538 	  )

   539 	 );

   541       while (!dfs.finished()) {

   542 	++dfs;

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

   544 	  if (dfs.isBNodeNewlyReached()) {

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

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

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

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

   551 	      free1.set

   552 		(w, std::min(free1[v], res_graph.resCap

   553 			     (typename ErasingResGW::Edge(dfs))));

   554 	    } else {

   555 	      free1.set

   556 		(w, res_graph.resCap

   557 		 (typename ErasingResGW::Edge(dfs)));

   558 	    }

   560 	    if (w==t) {

   561 	      __augment=true;

   562 	      _augment=true;

   563 	      break;

   564 	    }

   565 	  } else {

   566 	    erasing_res_graph.erase(dfs);

   567 	  }

   568 	}

   569       }

   571       if (__augment) {

   572 	typename ErasingResGW::Node

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

   574 	Num augment_value=free1[n];

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

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

   577 	  res_graph.augment(e, augment_value);

   578 	  n=erasing_res_graph.tail(e);

   579 	  if (res_graph.resCap(e)==0)

   580 	    erasing_res_graph.erase(e);

   581 	}

   582       }

   584     } //while (__augment)

   586     status=AFTER_AUGMENTING;

   587     return _augment;

   588   }

   591 } //namespace hugo

   593 #endif //HUGO_AUGMENTING_FLOW_H

author	marci
	Thu, 16 Sep 2004 13:54:46 +0000
changeset 865	2f3f87afb1d2
parent 862	732f2acb7239
child 888	cc3590763f7f
permissions	-rw-r--r--