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 Key;

   116       typedef bool Value;

   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->target(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->source(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.source(e);

   274 	Node w=res_graph.target(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.target(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.source(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.source(e);

   333 	Node w=res_graph.target(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.target(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.source(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       for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)

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

   387     }

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

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

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

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

   394     while ( !bfs.finished() ) {

   395       ResGWEdge e=bfs;

   396       if (e!=INVALID) {

   397 	if (bfs.isBNodeNewlyReached()) {

   398 	  dist.set(res_graph.target(e), dist[res_graph.source(e)]+1);

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

   400 					res_graph_to_F[res_graph.target(e)]);

   401 	  //original_edge.update();

   402 	  original_edge.set(f, e);

   403 	  //residual_capacity.update();

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

   405 	} else {

   406 	  if (dist[res_graph.target(e)]==(dist[res_graph.source(e)]+1)) {

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

   408 					  res_graph_to_F[res_graph.target(e)]);

   409 	    //original_edge.update();

   410 	    original_edge.set(f, e);

   411 	    //residual_capacity.update();

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

   413 	  }

   414 	}

   415       }

   416       ++bfs;

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

   419     bool __augment=true;

   421     while (__augment) {

   422       __augment=false;

   423       //computing blocking flow with dfs

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

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

   426       pred.set(sF, INVALID);

   427       //invalid iterators for sources

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

   431       dfs.pushAndSetReached(sF);

   432       while (!dfs.finished()) {

   433 	++dfs;

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

   435 	  if (dfs.isBNodeNewlyReached()) {

   436 	    typename MG::Node v=F.source(dfs);

   437 	    typename MG::Node w=F.target(dfs);

   438 	    pred.set(w, dfs);

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

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

   441 	    } else {

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

   443 	    }

   444 	    if (w==tF) {

   445 	      __augment=true;

   446 	      _augment=true;

   447 	      break;

   448 	    }

   450 	  } else {

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

   452 	  }

   453 	}

   454       }

   456       if (__augment) {

   457 	typename MG::Node n=tF;

   458 	Num augment_value=free[tF];

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

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

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

   462 	  n=F.source(e);

   463 	  if (residual_capacity[e]==augment_value)

   464 	    F.erase(e);

   465 	  else

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

   467 	}

   468       }

   470     }

   472     status=AFTER_AUGMENTING;

   473     return _augment;

   474   }

   476   /// Blocking flow augmentation without constructing the layered

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

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

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

   480   {

   481     bool _augment=false;

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

   484     typename ResGW::ResCap res_cap(res_graph);

   486     //Potential map, for distances from s

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

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

   489     Const1Map const_1_map(1);

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

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

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

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

   495     bfs.pushAndSetReached(s);

   497     //computing distances from s in the residual graph

   498     while ( !bfs.finished() ) {

   499       ResGWEdge e=bfs;

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

   501 	potential.set(res_graph.target(e), potential[res_graph.source(e)]+1);

   502       ++bfs;

   503     }

   505     //Subgraph containing the edges on some shortest paths

   506     //(i.e. tight edges)

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

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

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

   510       Const1Map> > FilterResGW;

   511     FilterResGW filter_res_graph(res_graph, true_map, tight_edge_filter);

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

   514     //met by the dfs

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

   516       first_out_edges(filter_res_graph);

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

   518       first_out_edges.set

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

   521     typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::

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

   523     ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);

   525     bool __augment=true;

   527     while (__augment) {

   529       __augment=false;

   530       //computing blocking flow with dfs

   531       DfsIterator< ErasingResGW,

   532 	typename ErasingResGW::template NodeMap<bool> >

   533 	dfs(erasing_res_graph);

   534       typename ErasingResGW::

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

   536       pred.set(s, INVALID);

   537       //invalid iterators for sources

   539       typename ErasingResGW::template NodeMap<Num>

   540 	free1(erasing_res_graph);

   542       dfs.pushAndSetReached

   543 	/// \bug lemon 0.2

   544 	(typename ErasingResGW::Node

   545 	 (typename FilterResGW::Node

   546 	  (typename ResGW::Node(s)

   547 	   )

   548 	  )

   549 	 );

   551       while (!dfs.finished()) {

   552 	++dfs;

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

   554 	  if (dfs.isBNodeNewlyReached()) {

   556 	    typename ErasingResGW::Node v=erasing_res_graph.source(dfs);

   557 	    typename ErasingResGW::Node w=erasing_res_graph.target(dfs);

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

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

   561 	      free1.set

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

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

   564 	    } else {

   565 	      free1.set

   566 		(w, res_cap

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

   568 	    }

   570 	    if (w==t) {

   571 	      __augment=true;

   572 	      _augment=true;

   573 	      break;

   574 	    }

   575 	  } else {

   576 	    erasing_res_graph.erase(dfs);

   577 	  }

   578 	}

   579       }

   581       if (__augment) {

   582 	typename ErasingResGW::Node

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

   584 	Num augment_value=free1[n];

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

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

   587 	  res_graph.augment(e, augment_value);

   588 	  n=erasing_res_graph.source(e);

   589 	  if (res_cap[e]==0)

   590 	    erasing_res_graph.erase(e);

   591 	}

   592       }

   594     } //while (__augment)

   596     status=AFTER_AUGMENTING;

   597     return _augment;

   598   }

   601 } //namespace lemon

   603 #endif //LEMON_AUGMENTING_FLOW_H

author	marci
	Wed, 17 Nov 2004 19:56:46 +0000
changeset 1004	b94037830dc8
parent 986	e997802b855c
permissions	-rw-r--r--