A bit better msg.
2 #ifndef LEMON_AUGMENTING_FLOW_H
3 #define LEMON_AUGMENTING_FLOW_H
8 #include <lemon/graph_wrapper.h>
11 #include <lemon/invalid.h>
12 #include <lemon/maps.h>
13 #include <demo/tight_edge_filter_map.h>
16 /// \brief Maximum flow algorithms.
20 using lemon::marci::BfsIterator;
21 using lemon::marci::DfsIterator;
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
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
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 {
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;
62 const CapMap* capacity;
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.
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
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.
99 AFTER_FAST_AUGMENTING,
100 AFTER_PRE_FLOW_PHASE_1,
101 AFTER_PRE_FLOW_PHASE_2
104 /// Don not needle this flag only if necessary.
106 int number_of_augmentations;
109 template<typename IntMap>
110 class TrickyReachedMap {
113 int* number_of_augmentations;
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) {
121 map->set(n, *number_of_augmentations);
123 map->set(n, *number_of_augmentations-1);
125 bool operator[](const Node& n) const {
126 return (*map)[n]==*number_of_augmentations;
130 AugmentingFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
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 {
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) {
170 case AFTER_PRE_FLOW_PHASE_2:
171 // std::cout << "AFTER_PRE_FLOW_PHASE_2" << std::endl;
174 // std::cout << "AFTER_NOTHING" << std::endl;
177 case AFTER_AUGMENTING:
178 // std::cout << "AFTER_AUGMENTING" << std::endl;
179 for(g->first(v); v!=INVALID; ++v) {
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) {
200 template<typename _CutMap>
201 void minMinCut(_CutMap& M) const {
202 std::queue<Node> queue;
207 while (!queue.empty()) {
208 Node w=queue.front();
212 for(g->first(e,w) ; e!=INVALID; ++e) {
214 if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
221 for(g->first(f,w) ; f!=INVALID; ++f) {
223 if (!M[v] && (*flow)[f] > 0 ) {
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;
239 Num flowValue() const {
241 for (InEdgeIt e(*g, t); e!=INVALID; ++e) a+=(*flow)[e];
242 for (OutEdgeIt e(*g, t); e!=INVALID; ++e) a-=(*flow)[e];
244 //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan
251 template <typename Graph, typename Num, typename CapMap, typename FlowMap>
252 bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()
254 ResGW res_graph(*g, *capacity, *flow);
255 typename ResGW::ResCap res_cap(res_graph);
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() ) {
272 if (e!=INVALID && bfs.isBNodeNewlyReached()) {
273 Node v=res_graph.source(e);
274 Node w=res_graph.target(e);
276 if (pred[v]!=INVALID) {
277 free.set(w, std::min(free[v], res_cap[e]));
279 free.set(w, res_cap[e]);
281 if (res_graph.target(e)==t) { _augment=true; break; }
285 } //end of searching augmenting path
289 Num augment_value=free[t];
290 while (pred[n]!=INVALID) {
292 res_graph.augment(e, augment_value);
293 n=res_graph.source(e);
297 status=AFTER_AUGMENTING;
301 template <typename Graph, typename Num, typename CapMap, typename FlowMap>
302 bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath2()
304 ResGW res_graph(*g, *capacity, *flow);
305 typename ResGW::ResCap res_cap(res_graph);
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;
313 ++number_of_augmentations;
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() ) {
331 if (e!=INVALID && bfs.isBNodeNewlyReached()) {
332 Node v=res_graph.source(e);
333 Node w=res_graph.target(e);
335 if (pred[v]!=INVALID) {
336 free.set(w, std::min(free[v], res_cap[e]));
338 free.set(w, res_cap[e]);
340 if (res_graph.target(e)==t) { _augment=true; break; }
344 } //end of searching augmenting path
348 Num augment_value=free[t];
349 while (pred[n]!=INVALID) {
351 res_graph.augment(e, augment_value);
352 n=res_graph.source(e);
356 status=AFTER_FAST_AUGMENTING;
361 template <typename Graph, typename Num, typename CapMap, typename FlowMap>
362 template<typename MutableGraph>
363 bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()
365 typedef MutableGraph MG;
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
382 typename ResGW::template NodeMap<typename MG::Node>
383 res_graph_to_F(res_graph);
385 for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)
386 res_graph_to_F.set(n, F.addNode());
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() ) {
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]);
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]);
417 } //computing distances from s in the residual graph
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()) {
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);
439 if (pred[v]!=INVALID) {
440 free.set(w, std::min(free[v], residual_capacity[dfs]));
442 free.set(w, residual_capacity[dfs]);
451 F.erase(typename MG::Edge(dfs));
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);
463 if (residual_capacity[e]==augment_value)
466 residual_capacity.set(e, residual_capacity[e]-augment_value);
472 status=AFTER_AUGMENTING;
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()
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() ) {
500 if (e!=INVALID && bfs.isBNodeNewlyReached())
501 potential.set(res_graph.target(e), potential[res_graph.source(e)]+1);
505 //Subgraph containing the edges on some shortest paths
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
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)
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);
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
544 (typename ErasingResGW::Node
545 (typename FilterResGW::Node
546 (typename ResGW::Node(s)
551 while (!dfs.finished()) {
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) {
562 (w, std::min(free1[v], res_cap
563 [typename ErasingResGW::Edge(dfs)]));
567 [typename ErasingResGW::Edge(dfs)]);
576 erasing_res_graph.erase(dfs);
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);
590 erasing_res_graph.erase(e);
594 } //while (__augment)
596 status=AFTER_AUGMENTING;
603 #endif //LEMON_AUGMENTING_FLOW_H