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

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

     2 #ifndef LEMON_AUGMENTING_FLOW_H

     3 #define LEMON_AUGMENTING_FLOW_H

     5 #include <vector>

     6 #include <iostream>

     8 #include <lemon/graph_wrapper.h>

     9 //#include <bfs_dfs.h>

    10 #include <bfs_mm.h>

    11 #include <lemon/invalid.h>

    12 #include <lemon/maps.h>

    13 #include <demo/tight_edge_filter_map.h>

    15 /// \file

    16 /// \brief Maximum flow algorithms.

    17 /// \ingroup galgs

    19 namespace lemon {

    20   using lemon::marci::BfsIterator;

    21   using lemon::marci::DfsIterator;

    23   /// \addtogroup galgs

    24   /// @{

    25   /// Class for augmenting path flow algorithms.

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

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

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

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

    31   /// constructor.

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

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

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

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

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

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

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

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

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

    42   /// cuts, resp.)

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

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

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

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

    47   ///\author Marton Makai

    48   template <typename Graph, typename Num,

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

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

    51   class AugmentingFlow {

    52   protected:

    53     typedef typename Graph::Node Node;

    54     typedef typename Graph::NodeIt NodeIt;

    55     typedef typename Graph::EdgeIt EdgeIt;

    56     typedef typename Graph::OutEdgeIt OutEdgeIt;

    57     typedef typename Graph::InEdgeIt InEdgeIt;

    59     const Graph* g;

    60     Node s;

    61     Node t;

    62     const CapMap* capacity;

    63     FlowMap* flow;

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

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

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

    67     typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;

    68     typedef typename ResGW::Edge ResGWEdge;

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

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

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

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

    74     //shows the levels of nodes.

    75     ReachedMap level;

    77   public:

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

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

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

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

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

    84     ///the \e target.

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

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

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

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

    89     enum FlowEnum{

    90       ZERO_FLOW,

    91       GEN_FLOW,

    92       PRE_FLOW,

    93       NO_FLOW

    94     };

    96     enum StatusEnum {

    97       AFTER_NOTHING,

    98       AFTER_AUGMENTING,

    99       AFTER_FAST_AUGMENTING,

   100       AFTER_PRE_FLOW_PHASE_1,

   101       AFTER_PRE_FLOW_PHASE_2

   102     };

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

   105     StatusEnum status;

   106     int number_of_augmentations;

   109     template<typename IntMap>

   110     class TrickyReachedMap {

   111     protected:

   112       IntMap* map;

   113       int* number_of_augmentations;

   114     public:

   115       typedef Node KeyType;

   116       typedef bool ValueType;

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

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

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

   120 	if (b)

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

   122 	else

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

   124       }

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

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

   127       }

   128     };

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

   131 		   FlowMap& _flow) :

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

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

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

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

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

   138     /// according to the Edmonds-Karp algorithm

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

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

   141     bool augmentOnShortestPath();

   142     bool augmentOnShortestPath2();

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

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

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

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

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

   149     template<typename MutableGraph> bool augmentOnBlockingFlow();

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

   152     /// residual graph is not constructed physically.

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

   154     bool augmentOnBlockingFlow2();

   156     template<typename _CutMap>

   157     void actMinCut(_CutMap& M) const {

   158       NodeIt v;

   159       switch (status) {

   160 	case AFTER_PRE_FLOW_PHASE_1:

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

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

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

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

   165 // 	  } else {

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

   167 // 	  }

   168 // 	}

   169 	break;

   170       case AFTER_PRE_FLOW_PHASE_2:

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

   172 	break;

   173       case AFTER_NOTHING:

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

   175 	minMinCut(M);

   176 	break;

   177       case AFTER_AUGMENTING:

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

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

   180 	  if (level[v]) {

   181 	    M.set(v, true);

   182 	  } else {

   183 	    M.set(v, false);

   184 	  }

   185 	}

   186 	break;

   187       case AFTER_FAST_AUGMENTING:

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

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

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

   191 	    M.set(v, true);

   192 	  } else {

   193 	    M.set(v, false);

   194 	  }

   195 	}

   196 	break;

   197       }

   198     }

   200     template<typename _CutMap>

   201     void minMinCut(_CutMap& M) const {

   202       std::queue<Node> queue;

   204       M.set(s,true);

   205       queue.push(s);

   207       while (!queue.empty()) {

   208         Node w=queue.front();

   209 	queue.pop();

   211 	OutEdgeIt e;

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

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

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

   215 	    queue.push(v);

   216 	    M.set(v, true);

   217 	  }

   218 	}

   220 	InEdgeIt f;

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

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

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

   224 	    queue.push(v);

   225 	    M.set(v, true);

   226 	  }

   227 	}

   228       }

   229     }

   231     template<typename _CutMap>

   232     void minMinCut2(_CutMap& M) const {

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

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

   235       bfs.pushAndSetReached(s);

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

   237     }

   239     Num flowValue() const {

   240       Num a=0;

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

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

   243       return a;

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

   245     }

   247   };

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

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

   253   {

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

   255     typename ResGW::ResCap res_cap(res_graph);

   257     bool _augment=false;

   259     //ReachedMap level(res_graph);

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

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

   262     bfs.pushAndSetReached(s);

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

   265     pred.set(s, INVALID);

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

   269     //searching for augmenting path

   270     while ( !bfs.finished() ) {

   271       ResGWEdge e=bfs;

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

   273 	Node v=res_graph.tail(e);

   274 	Node w=res_graph.head(e);

   275 	pred.set(w, e);

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

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

   278 	} else {

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

   280 	}

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

   282       }

   284       ++bfs;

   285     } //end of searching augmenting path

   287     if (_augment) {

   288       Node n=t;

   289       Num augment_value=free[t];

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

   291 	ResGWEdge e=pred[n];

   292 	res_graph.augment(e, augment_value);

   293 	n=res_graph.tail(e);

   294       }

   295     }

   297     status=AFTER_AUGMENTING;

   298     return _augment;

   299   }

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

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

   303   {

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

   305     typename ResGW::ResCap res_cap(res_graph);

   307     bool _augment=false;

   309     if (status!=AFTER_FAST_AUGMENTING) {

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

   311       number_of_augmentations=1;

   312     } else {

   313       ++number_of_augmentations;

   314     }

   315     TrickyReachedMap<ReachedMap>

   316       tricky_reached_map(level, number_of_augmentations);

   317     //ReachedMap level(res_graph);

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

   319     BfsIterator<ResGW, TrickyReachedMap<ReachedMap> >

   320       bfs(res_graph, tricky_reached_map);

   321     bfs.pushAndSetReached(s);

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

   324     pred.set(s, INVALID);

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

   328     //searching for augmenting path

   329     while ( !bfs.finished() ) {

   330       ResGWEdge e=bfs;

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

   332 	Node v=res_graph.tail(e);

   333 	Node w=res_graph.head(e);

   334 	pred.set(w, e);

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

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

   337 	} else {

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

   339 	}

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

   341       }

   343       ++bfs;

   344     } //end of searching augmenting path

   346     if (_augment) {

   347       Node n=t;

   348       Num augment_value=free[t];

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

   350 	ResGWEdge e=pred[n];

   351 	res_graph.augment(e, augment_value);

   352 	n=res_graph.tail(e);

   353       }

   354     }

   356     status=AFTER_FAST_AUGMENTING;

   357     return _augment;

   358   }

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

   362   template<typename MutableGraph>

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

   364   {

   365     typedef MutableGraph MG;

   366     bool _augment=false;

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

   369     typename ResGW::ResCap res_cap(res_graph);

   371     //bfs for distances on the residual graph

   372     //ReachedMap level(res_graph);

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

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

   375     bfs.pushAndSetReached(s);

   376     typename ResGW::template NodeMap<int>

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

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

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

   381     MG F;

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

   383       res_graph_to_F(res_graph);

   384     {

   385       typename ResGW::NodeIt n;

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

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

   388     }

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

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

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

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

   395     while ( !bfs.finished() ) {

   396       ResGWEdge e=bfs;

   397       if (e!=INVALID) {

   398 	if (bfs.isBNodeNewlyReached()) {

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

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

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

   402 	  //original_edge.update();

   403 	  original_edge.set(f, e);

   404 	  //residual_capacity.update();

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

   406 	} else {

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

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

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

   410 	    //original_edge.update();

   411 	    original_edge.set(f, e);

   412 	    //residual_capacity.update();

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

   414 	  }

   415 	}

   416       }

   417       ++bfs;

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

   420     bool __augment=true;

   422     while (__augment) {

   423       __augment=false;

   424       //computing blocking flow with dfs

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

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

   427       pred.set(sF, INVALID);

   428       //invalid iterators for sources

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

   432       dfs.pushAndSetReached(sF);

   433       while (!dfs.finished()) {

   434 	++dfs;

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

   436 	  if (dfs.isBNodeNewlyReached()) {

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

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

   439 	    pred.set(w, dfs);

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

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

   442 	    } else {

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

   444 	    }

   445 	    if (w==tF) {

   446 	      __augment=true;

   447 	      _augment=true;

   448 	      break;

   449 	    }

   451 	  } else {

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

   453 	  }

   454 	}

   455       }

   457       if (__augment) {

   458 	typename MG::Node n=tF;

   459 	Num augment_value=free[tF];

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

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

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

   463 	  n=F.tail(e);

   464 	  if (residual_capacity[e]==augment_value)

   465 	    F.erase(e);

   466 	  else

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

   468 	}

   469       }

   471     }

   473     status=AFTER_AUGMENTING;

   474     return _augment;

   475   }

   477   /// Blocking flow augmentation without constructing the layered

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

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

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

   481   {

   482     bool _augment=false;

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

   485     typename ResGW::ResCap res_cap(res_graph);

   487     //Potential map, for distances from s

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

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

   490     Const1Map const_1_map(1);

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

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

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

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

   496     bfs.pushAndSetReached(s);

   498     //computing distances from s in the residual graph

   499     while ( !bfs.finished() ) {

   500       ResGWEdge e=bfs;

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

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

   503       ++bfs;

   504     }

   506     //Subgraph containing the edges on some shortest paths

   507     //(i.e. tight edges)

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

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

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

   511       Const1Map> > FilterResGW;

   512     FilterResGW filter_res_graph(res_graph, true_map, tight_edge_filter);

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

   515     //met by the dfs

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

   517       first_out_edges(filter_res_graph);

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

   519       first_out_edges.set

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

   522     typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::

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

   524     ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);

   526     bool __augment=true;

   528     while (__augment) {

   530       __augment=false;

   531       //computing blocking flow with dfs

   532       DfsIterator< ErasingResGW,

   533 	typename ErasingResGW::template NodeMap<bool> >

   534 	dfs(erasing_res_graph);

   535       typename ErasingResGW::

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

   537       pred.set(s, INVALID);

   538       //invalid iterators for sources

   540       typename ErasingResGW::template NodeMap<Num>

   541 	free1(erasing_res_graph);

   543       dfs.pushAndSetReached

   544 	/// \bug lemon 0.2

   545 	(typename ErasingResGW::Node

   546 	 (typename FilterResGW::Node

   547 	  (typename ResGW::Node(s)

   548 	   )

   549 	  )

   550 	 );

   552       while (!dfs.finished()) {

   553 	++dfs;

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

   555 	  if (dfs.isBNodeNewlyReached()) {

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

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

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

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

   562 	      free1.set

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

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

   565 	    } else {

   566 	      free1.set

   567 		(w, res_cap

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

   569 	    }

   571 	    if (w==t) {

   572 	      __augment=true;

   573 	      _augment=true;

   574 	      break;

   575 	    }

   576 	  } else {

   577 	    erasing_res_graph.erase(dfs);

   578 	  }

   579 	}

   580       }

   582       if (__augment) {

   583 	typename ErasingResGW::Node

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

   585 	Num augment_value=free1[n];

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

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

   588 	  res_graph.augment(e, augment_value);

   589 	  n=erasing_res_graph.tail(e);

   590 	  if (res_cap[e]==0)

   591 	    erasing_res_graph.erase(e);

   592 	}

   593       }

   595     } //while (__augment)

   597     status=AFTER_AUGMENTING;

   598     return _augment;

   599   }

   602 } //namespace lemon

   604 #endif //LEMON_AUGMENTING_FLOW_H

author	alpar
	Mon, 08 Nov 2004 15:23:31 +0000
changeset 968	1a7593db0eaa
parent 921	818510fa3d99
child 970	09f9abe22df2
permissions	-rw-r--r--