marci@762: // -*- C++ -*- marci@762: #ifndef HUGO_AUGMENTING_FLOW_H marci@762: #define HUGO_AUGMENTING_FLOW_H marci@762: marci@762: #include marci@762: #include marci@762: marci@762: #include marci@762: #include marci@762: #include marci@762: #include marci@863: #include marci@762: marci@762: /// \file marci@762: /// \brief Maximum flow algorithms. marci@762: /// \ingroup galgs marci@762: marci@762: namespace hugo { marci@762: marci@762: /// \addtogroup galgs marci@762: /// @{ marci@862: /// Class for augmenting path flow algorithms. marci@762: marci@862: /// This class provides various algorithms for finding a flow of marci@862: /// maximum value in a directed graph. The \e source node, the \e marci@862: /// target node, the \e capacity of the edges and the \e starting \e marci@862: /// flow value of the edges should be passed to the algorithm through the marci@862: /// constructor. marci@862: // /// It is possible to change these quantities using the marci@862: // /// functions \ref resetSource, \ref resetTarget, \ref resetCap and marci@862: // /// \ref resetFlow. Before any subsequent runs of any algorithm of marci@862: // /// the class \ref resetFlow should be called. marci@762: marci@862: /// After running an algorithm of the class, the actual flow value marci@862: /// can be obtained by calling \ref flowValue(). The minimum marci@862: /// value cut can be written into a \c node map of \c bools by marci@862: /// calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes marci@862: /// the inclusionwise minimum and maximum of the minimum value marci@862: /// cuts, resp.) marci@762: ///\param Graph The directed graph type the algorithm runs on. marci@762: ///\param Num The number type of the capacities and the flow values. marci@762: ///\param CapMap The capacity map type. marci@762: ///\param FlowMap The flow map type. marci@862: ///\author Marton Makai marci@762: template , marci@762: typename FlowMap=typename Graph::template EdgeMap > marci@762: class AugmentingFlow { marci@762: protected: marci@762: typedef typename Graph::Node Node; marci@762: typedef typename Graph::NodeIt NodeIt; marci@762: typedef typename Graph::EdgeIt EdgeIt; marci@762: typedef typename Graph::OutEdgeIt OutEdgeIt; marci@762: typedef typename Graph::InEdgeIt InEdgeIt; marci@762: marci@762: const Graph* g; marci@762: Node s; marci@762: Node t; marci@762: const CapMap* capacity; marci@762: FlowMap* flow; marci@762: // int n; //the number of nodes of G marci@762: typedef ResGraphWrapper ResGW; marci@762: //typedef ExpResGraphWrapper ResGW; marci@762: typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt; marci@762: typedef typename ResGW::Edge ResGWEdge; marci@762: //typedef typename ResGW::template NodeMap ReachedMap; marci@762: typedef typename Graph::template NodeMap ReachedMap; marci@762: marci@762: //level works as a bool map in augmenting path algorithms and is marci@762: //used by bfs for storing reached information. In preflow, it marci@762: //shows the levels of nodes. marci@762: ReachedMap level; marci@762: marci@762: public: marci@762: ///Indicates the property of the starting flow. marci@762: marci@762: ///Indicates the property of the starting flow. The meanings are as follows: marci@762: ///- \c ZERO_FLOW: constant zero flow marci@762: ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to marci@762: ///the sum of the out-flows in every node except the \e source and marci@762: ///the \e target. marci@762: ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at marci@762: ///least the sum of the out-flows in every node except the \e source. marci@762: ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be marci@762: ///set to the constant zero flow in the beginning of the algorithm in this case. marci@762: enum FlowEnum{ marci@762: ZERO_FLOW, marci@762: GEN_FLOW, marci@762: PRE_FLOW, marci@762: NO_FLOW marci@762: }; marci@762: marci@762: enum StatusEnum { marci@762: AFTER_NOTHING, marci@762: AFTER_AUGMENTING, marci@762: AFTER_FAST_AUGMENTING, marci@762: AFTER_PRE_FLOW_PHASE_1, marci@762: AFTER_PRE_FLOW_PHASE_2 marci@762: }; marci@762: marci@762: /// Don not needle this flag only if necessary. marci@762: StatusEnum status; marci@762: int number_of_augmentations; marci@762: marci@762: marci@762: template marci@762: class TrickyReachedMap { marci@762: protected: marci@762: IntMap* map; marci@762: int* number_of_augmentations; marci@762: public: marci@762: TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) : marci@762: map(&_map), number_of_augmentations(&_number_of_augmentations) { } marci@762: void set(const Node& n, bool b) { marci@762: if (b) marci@762: map->set(n, *number_of_augmentations); marci@762: else marci@762: map->set(n, *number_of_augmentations-1); marci@762: } marci@762: bool operator[](const Node& n) const { marci@762: return (*map)[n]==*number_of_augmentations; marci@762: } marci@762: }; marci@762: marci@762: AugmentingFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, marci@762: FlowMap& _flow) : marci@762: g(&_G), s(_s), t(_t), capacity(&_capacity), marci@762: flow(&_flow), //n(_G.nodeNum()), marci@762: level(_G), //excess(_G,0), marci@762: status(AFTER_NOTHING), number_of_augmentations(0) { } marci@762: marci@762: /// Starting from a flow, this method searches for an augmenting path marci@762: /// according to the Edmonds-Karp algorithm marci@762: /// and augments the flow on if any. marci@762: /// The return value shows if the augmentation was succesful. marci@762: bool augmentOnShortestPath(); marci@762: bool augmentOnShortestPath2(); marci@762: marci@762: /// Starting from a flow, this method searches for an augmenting blocking marci@762: /// flow according to Dinits' algorithm and augments the flow on if any. marci@762: /// The blocking flow is computed in a physically constructed marci@762: /// residual graph of type \c Mutablegraph. marci@762: /// The return value show sif the augmentation was succesful. marci@762: template bool augmentOnBlockingFlow(); marci@762: marci@762: /// The same as \c augmentOnBlockingFlow but the marci@762: /// residual graph is not constructed physically. marci@762: /// The return value shows if the augmentation was succesful. marci@762: bool augmentOnBlockingFlow2(); marci@762: marci@762: template marci@762: void actMinCut(_CutMap& M) const { marci@762: NodeIt v; marci@762: switch (status) { marci@762: case AFTER_PRE_FLOW_PHASE_1: marci@762: // std::cout << "AFTER_PRE_FLOW_PHASE_1" << std::endl; marci@762: // for(g->first(v); g->valid(v); g->next(v)) { marci@762: // if (level[v] < n) { marci@762: // M.set(v, false); marci@762: // } else { marci@762: // M.set(v, true); marci@762: // } marci@762: // } marci@762: break; marci@762: case AFTER_PRE_FLOW_PHASE_2: marci@762: // std::cout << "AFTER_PRE_FLOW_PHASE_2" << std::endl; marci@762: break; marci@762: case AFTER_NOTHING: marci@762: // std::cout << "AFTER_NOTHING" << std::endl; marci@762: minMinCut(M); marci@762: break; marci@762: case AFTER_AUGMENTING: marci@762: // std::cout << "AFTER_AUGMENTING" << std::endl; marci@775: for(g->first(v); v!=INVALID; ++v) { marci@762: if (level[v]) { marci@762: M.set(v, true); marci@762: } else { marci@762: M.set(v, false); marci@762: } marci@762: } marci@762: break; marci@762: case AFTER_FAST_AUGMENTING: marci@762: // std::cout << "AFTER_FAST_AUGMENTING" << std::endl; marci@775: for(g->first(v); v!=INVALID; ++v) { marci@762: if (level[v]==number_of_augmentations) { marci@762: M.set(v, true); marci@762: } else { marci@762: M.set(v, false); marci@762: } marci@762: } marci@762: break; marci@762: } marci@762: } marci@762: marci@762: template marci@762: void minMinCut(_CutMap& M) const { marci@762: std::queue queue; marci@762: marci@762: M.set(s,true); marci@762: queue.push(s); marci@762: marci@762: while (!queue.empty()) { marci@762: Node w=queue.front(); marci@762: queue.pop(); marci@762: marci@762: OutEdgeIt e; marci@775: for(g->first(e,w) ; e!=INVALID; ++e) { marci@762: Node v=g->head(e); marci@762: if (!M[v] && (*flow)[e] < (*capacity)[e] ) { marci@762: queue.push(v); marci@762: M.set(v, true); marci@762: } marci@762: } marci@762: marci@762: InEdgeIt f; marci@775: for(g->first(f,w) ; f!=INVALID; ++f) { marci@762: Node v=g->tail(f); marci@762: if (!M[v] && (*flow)[f] > 0 ) { marci@762: queue.push(v); marci@762: M.set(v, true); marci@762: } marci@762: } marci@762: } marci@762: } marci@762: marci@762: template marci@762: void minMinCut2(_CutMap& M) const { marci@762: ResGW res_graph(*g, *capacity, *flow); marci@762: BfsIterator bfs(res_graph, M); marci@762: bfs.pushAndSetReached(s); marci@762: while (!bfs.finished()) ++bfs; marci@762: } marci@762: marci@762: Num flowValue() const { marci@762: Num a=0; marci@777: for (InEdgeIt e(*g, t); e!=INVALID; ++e) a+=(*flow)[e]; marci@777: for (OutEdgeIt e(*g, t); e!=INVALID; ++e) a-=(*flow)[e]; marci@762: return a; marci@762: //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan marci@762: } marci@762: marci@762: }; marci@762: marci@762: marci@762: marci@762: template marci@762: bool AugmentingFlow::augmentOnShortestPath() marci@762: { marci@762: ResGW res_graph(*g, *capacity, *flow); marci@762: bool _augment=false; marci@762: marci@762: //ReachedMap level(res_graph); marci@777: for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0); marci@762: BfsIterator bfs(res_graph, level); marci@762: bfs.pushAndSetReached(s); marci@762: marci@762: typename ResGW::template NodeMap pred(res_graph); marci@762: pred.set(s, INVALID); marci@762: marci@762: typename ResGW::template NodeMap free(res_graph); marci@762: marci@762: //searching for augmenting path marci@762: while ( !bfs.finished() ) { marci@777: ResGWEdge e=bfs; marci@775: if (e!=INVALID && bfs.isBNodeNewlyReached()) { marci@762: Node v=res_graph.tail(e); marci@762: Node w=res_graph.head(e); marci@762: pred.set(w, e); marci@775: if (pred[v]!=INVALID) { marci@762: free.set(w, std::min(free[v], res_graph.resCap(e))); marci@762: } else { marci@762: free.set(w, res_graph.resCap(e)); marci@762: } marci@762: if (res_graph.head(e)==t) { _augment=true; break; } marci@762: } marci@762: marci@762: ++bfs; marci@762: } //end of searching augmenting path marci@762: marci@762: if (_augment) { marci@762: Node n=t; marci@762: Num augment_value=free[t]; marci@775: while (pred[n]!=INVALID) { marci@762: ResGWEdge e=pred[n]; marci@762: res_graph.augment(e, augment_value); marci@762: n=res_graph.tail(e); marci@762: } marci@762: } marci@762: marci@762: status=AFTER_AUGMENTING; marci@762: return _augment; marci@762: } marci@762: marci@762: template marci@762: bool AugmentingFlow::augmentOnShortestPath2() marci@762: { marci@762: ResGW res_graph(*g, *capacity, *flow); marci@762: bool _augment=false; marci@762: marci@762: if (status!=AFTER_FAST_AUGMENTING) { marci@777: for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0); marci@762: number_of_augmentations=1; marci@762: } else { marci@762: ++number_of_augmentations; marci@762: } marci@762: TrickyReachedMap marci@762: tricky_reached_map(level, number_of_augmentations); marci@762: //ReachedMap level(res_graph); marci@762: // FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); marci@762: BfsIterator > marci@762: bfs(res_graph, tricky_reached_map); marci@762: bfs.pushAndSetReached(s); marci@762: marci@762: typename ResGW::template NodeMap pred(res_graph); marci@762: pred.set(s, INVALID); marci@762: marci@762: typename ResGW::template NodeMap free(res_graph); marci@762: marci@762: //searching for augmenting path marci@762: while ( !bfs.finished() ) { marci@777: ResGWEdge e=bfs; marci@775: if (e!=INVALID && bfs.isBNodeNewlyReached()) { marci@762: Node v=res_graph.tail(e); marci@762: Node w=res_graph.head(e); marci@762: pred.set(w, e); marci@775: if (pred[v]!=INVALID) { marci@762: free.set(w, std::min(free[v], res_graph.resCap(e))); marci@762: } else { marci@762: free.set(w, res_graph.resCap(e)); marci@762: } marci@762: if (res_graph.head(e)==t) { _augment=true; break; } marci@762: } marci@762: marci@762: ++bfs; marci@762: } //end of searching augmenting path marci@762: marci@762: if (_augment) { marci@762: Node n=t; marci@762: Num augment_value=free[t]; marci@775: while (pred[n]!=INVALID) { marci@762: ResGWEdge e=pred[n]; marci@762: res_graph.augment(e, augment_value); marci@762: n=res_graph.tail(e); marci@762: } marci@762: } marci@762: marci@762: status=AFTER_FAST_AUGMENTING; marci@762: return _augment; marci@762: } marci@762: marci@762: marci@762: template marci@762: template marci@762: bool AugmentingFlow::augmentOnBlockingFlow() marci@762: { marci@762: typedef MutableGraph MG; marci@762: bool _augment=false; marci@762: marci@762: ResGW res_graph(*g, *capacity, *flow); marci@762: marci@762: //bfs for distances on the residual graph marci@762: //ReachedMap level(res_graph); marci@777: for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0); marci@762: BfsIterator bfs(res_graph, level); marci@762: bfs.pushAndSetReached(s); marci@762: typename ResGW::template NodeMap marci@762: dist(res_graph); //filled up with 0's marci@762: marci@762: //F will contain the physical copy of the residual graph marci@762: //with the set of edges which are on shortest paths marci@762: MG F; marci@762: typename ResGW::template NodeMap marci@762: res_graph_to_F(res_graph); marci@762: { marci@762: typename ResGW::NodeIt n; marci@862: for(res_graph.first(n); n!=INVALID; ++n) marci@762: res_graph_to_F.set(n, F.addNode()); marci@762: } marci@762: marci@762: typename MG::Node sF=res_graph_to_F[s]; marci@762: typename MG::Node tF=res_graph_to_F[t]; marci@762: typename MG::template EdgeMap original_edge(F); marci@762: typename MG::template EdgeMap residual_capacity(F); marci@762: marci@762: while ( !bfs.finished() ) { marci@777: ResGWEdge e=bfs; marci@775: if (e!=INVALID) { marci@762: if (bfs.isBNodeNewlyReached()) { marci@762: dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1); marci@762: typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], marci@762: res_graph_to_F[res_graph.head(e)]); marci@854: //original_edge.update(); marci@762: original_edge.set(f, e); marci@854: //residual_capacity.update(); marci@762: residual_capacity.set(f, res_graph.resCap(e)); marci@762: } else { marci@762: if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) { marci@762: typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], marci@762: res_graph_to_F[res_graph.head(e)]); marci@854: //original_edge.update(); marci@762: original_edge.set(f, e); marci@854: //residual_capacity.update(); marci@762: residual_capacity.set(f, res_graph.resCap(e)); marci@762: } marci@762: } marci@762: } marci@762: ++bfs; marci@762: } //computing distances from s in the residual graph marci@762: marci@762: bool __augment=true; marci@762: marci@762: while (__augment) { marci@762: __augment=false; marci@762: //computing blocking flow with dfs marci@762: DfsIterator< MG, typename MG::template NodeMap > dfs(F); marci@762: typename MG::template NodeMap pred(F); marci@762: pred.set(sF, INVALID); marci@762: //invalid iterators for sources marci@762: marci@762: typename MG::template NodeMap free(F); marci@762: marci@762: dfs.pushAndSetReached(sF); marci@762: while (!dfs.finished()) { marci@762: ++dfs; marci@854: if (typename MG::Edge(dfs)!=INVALID) { marci@762: if (dfs.isBNodeNewlyReached()) { marci@777: typename MG::Node v=F.tail(dfs); marci@777: typename MG::Node w=F.head(dfs); marci@762: pred.set(w, dfs); marci@775: if (pred[v]!=INVALID) { marci@762: free.set(w, std::min(free[v], residual_capacity[dfs])); marci@762: } else { marci@762: free.set(w, residual_capacity[dfs]); marci@762: } marci@762: if (w==tF) { marci@762: __augment=true; marci@762: _augment=true; marci@762: break; marci@762: } marci@762: marci@762: } else { marci@854: F.erase(typename MG::Edge(dfs)); marci@762: } marci@762: } marci@762: } marci@762: marci@762: if (__augment) { marci@762: typename MG::Node n=tF; marci@762: Num augment_value=free[tF]; marci@775: while (pred[n]!=INVALID) { marci@762: typename MG::Edge e=pred[n]; marci@762: res_graph.augment(original_edge[e], augment_value); marci@762: n=F.tail(e); marci@762: if (residual_capacity[e]==augment_value) marci@762: F.erase(e); marci@762: else marci@762: residual_capacity.set(e, residual_capacity[e]-augment_value); marci@762: } marci@762: } marci@762: marci@762: } marci@762: marci@762: status=AFTER_AUGMENTING; marci@762: return _augment; marci@762: } marci@762: marci@862: /// Blocking flow augmentation without constructing the layered marci@862: /// graph physically in which the blocking flow is computed. marci@762: template marci@762: bool AugmentingFlow::augmentOnBlockingFlow2() marci@762: { marci@762: bool _augment=false; marci@762: marci@762: ResGW res_graph(*g, *capacity, *flow); marci@762: marci@862: //Potential map, for distances from s marci@862: typename ResGW::template NodeMap potential(res_graph, 0); marci@862: typedef ConstMap Const1Map; marci@862: Const1Map const_1_map(1); marci@862: TightEdgeFilterMap, marci@862: Const1Map> tight_edge_filter(res_graph, potential, const_1_map); marci@862: marci@777: for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0); marci@762: BfsIterator bfs(res_graph, level); marci@862: bfs.pushAndSetReached(s); marci@762: marci@862: //computing distances from s in the residual graph marci@762: while ( !bfs.finished() ) { marci@777: ResGWEdge e=bfs; marci@862: if (e!=INVALID && bfs.isBNodeNewlyReached()) marci@862: potential.set(res_graph.head(e), potential[res_graph.tail(e)]+1); marci@762: ++bfs; marci@862: } marci@762: marci@862: //Subgraph containing the edges on some shortest paths marci@862: //(i.e. tight edges) marci@762: ConstMap true_map(true); marci@762: typedef SubGraphWrapper, marci@862: TightEdgeFilterMap, marci@862: Const1Map> > FilterResGW; marci@862: FilterResGW filter_res_graph(res_graph, true_map, tight_edge_filter); marci@762: marci@762: //Subgraph, which is able to delete edges which are already marci@762: //met by the dfs marci@777: typename FilterResGW::template NodeMap marci@762: first_out_edges(filter_res_graph); marci@862: for (typename FilterResGW::NodeIt v(filter_res_graph); v!=INVALID; ++v) marci@862: first_out_edges.set marci@862: (v, typename FilterResGW::OutEdgeIt(filter_res_graph, v)); marci@862: marci@762: typedef ErasingFirstGraphWrapper > ErasingResGW; marci@762: ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges); marci@762: marci@762: bool __augment=true; marci@762: marci@762: while (__augment) { marci@762: marci@762: __augment=false; marci@762: //computing blocking flow with dfs marci@762: DfsIterator< ErasingResGW, marci@762: typename ErasingResGW::template NodeMap > marci@762: dfs(erasing_res_graph); marci@762: typename ErasingResGW:: marci@777: template NodeMap pred(erasing_res_graph); marci@762: pred.set(s, INVALID); marci@762: //invalid iterators for sources marci@762: marci@762: typename ErasingResGW::template NodeMap marci@762: free1(erasing_res_graph); marci@762: marci@762: dfs.pushAndSetReached marci@777: /// \bug hugo 0.2 marci@762: (typename ErasingResGW::Node marci@762: (typename FilterResGW::Node marci@762: (typename ResGW::Node(s) marci@762: ) marci@762: ) marci@762: ); marci@777: marci@762: while (!dfs.finished()) { marci@762: ++dfs; marci@862: if (typename ErasingResGW::Edge(dfs)!=INVALID) { marci@862: if (dfs.isBNodeNewlyReached()) { marci@862: marci@862: typename ErasingResGW::Node v=erasing_res_graph.tail(dfs); marci@862: typename ErasingResGW::Node w=erasing_res_graph.head(dfs); marci@762: marci@862: pred.set(w, typename ErasingResGW::Edge(dfs)); marci@862: if (pred[v]!=INVALID) { marci@862: free1.set marci@862: (w, std::min(free1[v], res_graph.resCap marci@862: (typename ErasingResGW::Edge(dfs)))); marci@862: } else { marci@862: free1.set marci@862: (w, res_graph.resCap marci@862: (typename ErasingResGW::Edge(dfs))); marci@862: } marci@762: marci@862: if (w==t) { marci@862: __augment=true; marci@862: _augment=true; marci@862: break; marci@762: } marci@862: } else { marci@862: erasing_res_graph.erase(dfs); marci@762: } marci@862: } marci@762: } marci@762: marci@762: if (__augment) { marci@762: typename ErasingResGW::Node marci@762: n=typename FilterResGW::Node(typename ResGW::Node(t)); marci@762: Num augment_value=free1[n]; marci@777: while (pred[n]!=INVALID) { marci@777: typename ErasingResGW::Edge e=pred[n]; marci@762: res_graph.augment(e, augment_value); marci@762: n=erasing_res_graph.tail(e); marci@762: if (res_graph.resCap(e)==0) marci@762: erasing_res_graph.erase(e); marci@762: } marci@762: } marci@762: marci@762: } //while (__augment) marci@762: marci@762: status=AFTER_AUGMENTING; marci@762: return _augment; marci@762: } marci@762: marci@762: marci@762: } //namespace hugo marci@762: marci@762: #endif //HUGO_AUGMENTING_FLOW_H marci@762: marci@762: