marci@478: // -*- C++ -*- marci@480: #ifndef HUGO_MAX_FLOW_H marci@480: #define HUGO_MAX_FLOW_H marci@478: marci@478: #include marci@478: #include marci@478: #include marci@478: marci@557: #include marci@602: #include marci@555: #include marci@555: #include marci@640: #include marci@478: marci@488: /// \file jacint@631: /// \brief Maximum flow algorithms. marci@615: /// \ingroup galgs marci@478: marci@478: namespace hugo { marci@478: jacint@631: /// \addtogroup galgs jacint@631: /// @{ jacint@631: ///Maximum flow algorithms class. marci@488: jacint@631: ///This class provides various algorithms for finding a flow of jacint@631: ///maximum value in a directed graph. The \e source node, the \e jacint@631: ///target node, the \e capacity of the edges and the \e starting \e marci@647: ///flow value of the edges should be passed to the algorithm through the jacint@631: ///constructor. It is possible to change these quantities using the jacint@631: ///functions \ref resetSource, \ref resetTarget, \ref resetCap and jacint@631: ///\ref resetFlow. Before any subsequent runs of any algorithm of marci@647: ///the class \ref resetFlow should be called. marci@647: marci@647: ///After running an algorithm of the class, the actual flow value marci@647: ///can be obtained by calling \ref flowValue(). The minimum jacint@631: ///value cut can be written into a \c node map of \c bools by jacint@631: ///calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes jacint@631: ///the inclusionwise minimum and maximum of the minimum value jacint@631: ///cuts, resp.) marci@632: ///\param Graph The directed graph type the algorithm runs on. jacint@631: ///\param Num The number type of the capacities and the flow values. marci@647: ///\param CapMap The capacity map type. marci@647: ///\param FlowMap The flow map type. jacint@631: ///\author Marton Makai, Jacint Szabo marci@615: template , marci@478: typename FlowMap=typename Graph::template EdgeMap > marci@478: class MaxFlow { marci@615: protected: marci@478: typedef typename Graph::Node Node; marci@478: typedef typename Graph::NodeIt NodeIt; jacint@631: typedef typename Graph::EdgeIt EdgeIt; marci@478: typedef typename Graph::OutEdgeIt OutEdgeIt; marci@478: typedef typename Graph::InEdgeIt InEdgeIt; marci@478: marci@478: typedef typename std::vector > VecStack; marci@478: typedef typename Graph::template NodeMap NNMap; marci@478: typedef typename std::vector VecNode; marci@478: marci@478: const Graph* g; marci@478: Node s; marci@478: Node t; marci@615: const CapMap* capacity; marci@478: FlowMap* flow; marci@478: int n; //the number of nodes of G marci@653: typedef ResGraphWrapper ResGW; marci@653: //typedef ExpResGraphWrapper ResGW; marci@478: typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt; marci@478: typedef typename ResGW::Edge ResGWEdge; marci@478: //typedef typename ResGW::template NodeMap ReachedMap; marci@478: typedef typename Graph::template NodeMap ReachedMap; jacint@631: jacint@631: jacint@631: //level works as a bool map in augmenting path algorithms and is jacint@631: //used by bfs for storing reached information. In preflow, it jacint@631: //shows the levels of nodes. marci@478: ReachedMap level; jacint@631: jacint@631: //excess is needed only in preflow marci@615: typename Graph::template NodeMap excess; jacint@631: jacint@631: //fixme jacint@631: // protected: marci@602: // MaxFlow() { } marci@615: // void set(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, marci@615: // FlowMap& _flow) marci@602: // { marci@615: // g=&_G; marci@615: // s=_s; marci@615: // t=_t; marci@602: // capacity=&_capacity; marci@602: // flow=&_flow; marci@602: // n=_G.nodeNum; marci@615: // level.set (_G); //kellene vmi ilyesmi fv marci@602: // excess(_G,0); //itt is marci@602: // } marci@478: marci@615: // constants used for heuristics marci@615: static const int H0=20; marci@615: static const int H1=1; marci@615: marci@478: public: marci@615: jacint@631: ///Indicates the property of the starting flow. jacint@631: jacint@631: ///Indicates the property of the starting flow. The meanings are as follows: jacint@631: ///- \c ZERO_FLOW: constant zero flow jacint@631: ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to jacint@631: ///the sum of the out-flows in every node except the \e source and jacint@631: ///the \e target. jacint@631: ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at jacint@631: ///least the sum of the out-flows in every node except the \e source. jacint@631: ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be jacint@631: ///set to the constant zero flow in the beginning of the algorithm in this case. marci@647: enum FlowEnum{ marci@615: ZERO_FLOW, marci@615: GEN_FLOW, marci@615: PRE_FLOW, marci@615: NO_FLOW marci@478: }; marci@478: marci@647: enum StatusEnum { marci@647: AFTER_NOTHING, marci@647: AFTER_AUGMENTING, marci@656: AFTER_FAST_AUGMENTING, marci@647: AFTER_PRE_FLOW_PHASE_1, marci@647: AFTER_PRE_FLOW_PHASE_2 marci@647: }; marci@647: marci@647: /// Don not needle this flag only if necessary. marci@647: StatusEnum status; marci@647: int number_of_augmentations; marci@647: marci@647: marci@647: template marci@647: class TrickyReachedMap { marci@647: protected: marci@647: IntMap* map; marci@647: int* number_of_augmentations; marci@647: public: marci@647: TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) : marci@647: map(&_map), number_of_augmentations(&_number_of_augmentations) { } marci@647: void set(const Node& n, bool b) { marci@650: if (b) marci@650: map->set(n, *number_of_augmentations); marci@650: else marci@650: map->set(n, *number_of_augmentations-1); marci@647: } marci@647: bool operator[](const Node& n) const { marci@647: return (*map)[n]==*number_of_augmentations; marci@647: } marci@647: }; marci@647: alpar@709: ///Constructor alpar@709: alpar@709: ///\todo Document, please. alpar@709: /// marci@615: MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, marci@478: FlowMap& _flow) : marci@615: g(&_G), s(_s), t(_t), capacity(&_capacity), marci@647: flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0), marci@647: status(AFTER_NOTHING), number_of_augmentations(0) { } marci@478: jacint@631: ///Runs a maximum flow algorithm. jacint@631: jacint@631: ///Runs a preflow algorithm, which is the fastest maximum flow jacint@631: ///algorithm up-to-date. The default for \c fe is ZERO_FLOW. jacint@631: ///\pre The starting flow must be jacint@631: /// - a constant zero flow if \c fe is \c ZERO_FLOW, jacint@631: /// - an arbitary flow if \c fe is \c GEN_FLOW, jacint@631: /// - an arbitary preflow if \c fe is \c PRE_FLOW, jacint@631: /// - any map if \c fe is NO_FLOW. marci@647: void run(FlowEnum fe=ZERO_FLOW) { marci@615: preflow(fe); marci@478: } marci@615: marci@647: jacint@631: ///Runs a preflow algorithm. jacint@631: jacint@631: ///Runs a preflow algorithm. The preflow algorithms provide the jacint@631: ///fastest way to compute a maximum flow in a directed graph. jacint@631: ///\pre The starting flow must be jacint@631: /// - a constant zero flow if \c fe is \c ZERO_FLOW, jacint@631: /// - an arbitary flow if \c fe is \c GEN_FLOW, jacint@631: /// - an arbitary preflow if \c fe is \c PRE_FLOW, jacint@631: /// - any map if \c fe is NO_FLOW. alpar@709: /// alpar@709: ///\todo NO_FLOW should be the default flow. marci@647: void preflow(FlowEnum fe) { jacint@631: preflowPhase1(fe); jacint@631: preflowPhase2(); marci@478: } jacint@631: // Heuristics: jacint@631: // 2 phase jacint@631: // gap jacint@631: // list 'level_list' on the nodes on level i implemented by hand jacint@631: // stack 'active' on the active nodes on level i jacint@631: // runs heuristic 'highest label' for H1*n relabels jacint@631: // runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label' jacint@631: // Parameters H0 and H1 are initialized to 20 and 1. marci@478: jacint@631: ///Runs the first phase of the preflow algorithm. marci@478: jacint@631: ///The preflow algorithm consists of two phases, this method runs the jacint@631: ///first phase. After the first phase the maximum flow value and a jacint@631: ///minimum value cut can already be computed, though a maximum flow jacint@631: ///is net yet obtained. So after calling this method \ref flowValue jacint@631: ///and \ref actMinCut gives proper results. jacint@631: ///\warning: \ref minCut, \ref minMinCut and \ref maxMinCut do not jacint@631: ///give minimum value cuts unless calling \ref preflowPhase2. jacint@631: ///\pre The starting flow must be jacint@631: /// - a constant zero flow if \c fe is \c ZERO_FLOW, jacint@631: /// - an arbitary flow if \c fe is \c GEN_FLOW, jacint@631: /// - an arbitary preflow if \c fe is \c PRE_FLOW, jacint@631: /// - any map if \c fe is NO_FLOW. marci@647: void preflowPhase1(FlowEnum fe); jacint@631: jacint@631: ///Runs the second phase of the preflow algorithm. jacint@631: jacint@631: ///The preflow algorithm consists of two phases, this method runs jacint@631: ///the second phase. After calling \ref preflowPhase1 and then jacint@631: ///\ref preflowPhase2 the methods \ref flowValue, \ref minCut, jacint@631: ///\ref minMinCut and \ref maxMinCut give proper results. jacint@631: ///\pre \ref preflowPhase1 must be called before. jacint@631: void preflowPhase2(); marci@478: marci@615: /// Starting from a flow, this method searches for an augmenting path marci@615: /// according to the Edmonds-Karp algorithm marci@615: /// and augments the flow on if any. marci@487: /// The return value shows if the augmentation was succesful. marci@478: bool augmentOnShortestPath(); marci@647: bool augmentOnShortestPath2(); marci@478: marci@615: /// Starting from a flow, this method searches for an augmenting blocking marci@615: /// flow according to Dinits' algorithm and augments the flow on if any. marci@615: /// The blocking flow is computed in a physically constructed marci@485: /// residual graph of type \c Mutablegraph. marci@487: /// The return value show sif the augmentation was succesful. marci@478: template bool augmentOnBlockingFlow(); marci@478: marci@615: /// The same as \c augmentOnBlockingFlow but the marci@485: /// residual graph is not constructed physically. marci@487: /// The return value shows if the augmentation was succesful. marci@478: bool augmentOnBlockingFlow2(); marci@478: jacint@631: /// Returns the maximum value of a flow. jacint@631: jacint@631: /// Returns the maximum value of a flow, by counting the jacint@631: /// over-flow of the target node \ref t. jacint@631: /// It can be called already after running \ref preflowPhase1. marci@647: Num flowValue() const { marci@478: Num a=0; marci@615: FOR_EACH_INC_LOC(InEdgeIt, e, *g, t) a+=(*flow)[e]; marci@615: FOR_EACH_INC_LOC(OutEdgeIt, e, *g, t) a-=(*flow)[e]; marci@478: return a; jacint@631: //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan marci@478: } marci@478: jacint@631: ///Returns a minimum value cut after calling \ref preflowPhase1. jacint@631: jacint@631: ///After the first phase of the preflow algorithm the maximum flow jacint@631: ///value and a minimum value cut can already be computed. This jacint@631: ///method can be called after running \ref preflowPhase1 for jacint@631: ///obtaining a minimum value cut. jacint@631: /// \warning Gives proper result only right after calling \ref jacint@631: /// preflowPhase1. marci@615: /// \todo We have to make some status variable which shows the marci@615: /// actual state marci@615: /// of the class. This enables us to determine which methods are valid marci@485: /// for MinCut computation marci@478: template marci@647: void actMinCut(_CutMap& M) const { marci@478: NodeIt v; marci@647: switch (status) { marci@656: case AFTER_PRE_FLOW_PHASE_1: marci@647: for(g->first(v); g->valid(v); g->next(v)) { marci@647: if (level[v] < n) { marci@647: M.set(v, false); marci@647: } else { marci@647: M.set(v, true); marci@647: } marci@485: } marci@647: break; marci@656: case AFTER_PRE_FLOW_PHASE_2: marci@656: case AFTER_NOTHING: marci@647: minMinCut(M); marci@647: break; marci@656: case AFTER_AUGMENTING: marci@647: for(g->first(v); g->valid(v); g->next(v)) { marci@647: if (level[v]) { marci@647: M.set(v, true); marci@647: } else { marci@647: M.set(v, false); marci@647: } marci@647: } marci@647: break; marci@656: case AFTER_FAST_AUGMENTING: marci@656: for(g->first(v); g->valid(v); g->next(v)) { marci@656: if (level[v]==number_of_augmentations) { marci@656: M.set(v, true); marci@656: } else { marci@656: M.set(v, false); marci@656: } marci@656: } marci@656: break; marci@478: } marci@478: } marci@478: jacint@631: ///Returns the inclusionwise minimum of the minimum value cuts. jacint@631: jacint@631: ///Sets \c M to the characteristic vector of the minimum value cut jacint@631: ///which is inclusionwise minimum. It is computed by processing jacint@631: ///a bfs from the source node \c s in the residual graph. jacint@631: ///\pre M should be a node map of bools initialized to false. jacint@631: ///\pre \c flow must be a maximum flow. marci@478: template marci@647: void minMinCut(_CutMap& M) const { marci@478: std::queue queue; marci@615: marci@615: M.set(s,true); marci@478: queue.push(s); marci@478: marci@478: while (!queue.empty()) { marci@478: Node w=queue.front(); marci@478: queue.pop(); marci@478: marci@478: OutEdgeIt e; marci@478: for(g->first(e,w) ; g->valid(e); g->next(e)) { marci@478: Node v=g->head(e); marci@478: if (!M[v] && (*flow)[e] < (*capacity)[e] ) { marci@478: queue.push(v); marci@478: M.set(v, true); marci@478: } marci@615: } marci@478: marci@478: InEdgeIt f; marci@478: for(g->first(f,w) ; g->valid(f); g->next(f)) { marci@478: Node v=g->tail(f); marci@478: if (!M[v] && (*flow)[f] > 0 ) { marci@478: queue.push(v); marci@478: M.set(v, true); marci@478: } marci@615: } marci@478: } marci@478: } marci@478: jacint@631: ///Returns the inclusionwise maximum of the minimum value cuts. marci@478: jacint@631: ///Sets \c M to the characteristic vector of the minimum value cut jacint@631: ///which is inclusionwise maximum. It is computed by processing a jacint@631: ///backward bfs from the target node \c t in the residual graph. jacint@631: ///\pre M should be a node map of bools initialized to false. jacint@631: ///\pre \c flow must be a maximum flow. marci@478: template marci@647: void maxMinCut(_CutMap& M) const { marci@478: marci@478: NodeIt v; marci@478: for(g->first(v) ; g->valid(v); g->next(v)) { marci@478: M.set(v, true); marci@478: } marci@478: marci@478: std::queue queue; marci@615: marci@615: M.set(t,false); marci@478: queue.push(t); marci@478: marci@478: while (!queue.empty()) { marci@478: Node w=queue.front(); marci@478: queue.pop(); marci@478: marci@478: InEdgeIt e; marci@478: for(g->first(e,w) ; g->valid(e); g->next(e)) { marci@478: Node v=g->tail(e); marci@478: if (M[v] && (*flow)[e] < (*capacity)[e] ) { marci@478: queue.push(v); marci@478: M.set(v, false); marci@478: } marci@478: } marci@615: marci@478: OutEdgeIt f; marci@478: for(g->first(f,w) ; g->valid(f); g->next(f)) { marci@478: Node v=g->head(f); marci@478: if (M[v] && (*flow)[f] > 0 ) { marci@478: queue.push(v); marci@478: M.set(v, false); marci@478: } marci@478: } marci@478: } marci@478: } marci@478: jacint@631: ///Returns a minimum value cut. marci@478: jacint@631: ///Sets \c M to the characteristic vector of a minimum value cut. jacint@631: ///\pre M should be a node map of bools initialized to false. jacint@631: ///\pre \c flow must be a maximum flow. marci@478: template marci@647: void minCut(CutMap& M) const { minMinCut(M); } marci@478: jacint@631: ///Resets the source node to \c _s. jacint@631: jacint@631: ///Resets the source node to \c _s. jacint@631: /// marci@647: void resetSource(Node _s) { s=_s; status=AFTER_NOTHING; } jacint@631: jacint@631: ///Resets the target node to \c _t. jacint@631: jacint@631: ///Resets the target node to \c _t. marci@487: /// marci@647: void resetTarget(Node _t) { t=_t; status=AFTER_NOTHING; } marci@615: jacint@631: /// Resets the edge map of the capacities to _cap. jacint@631: jacint@631: /// Resets the edge map of the capacities to _cap. jacint@631: /// marci@647: void resetCap(const CapMap& _cap) { capacity=&_cap; status=AFTER_NOTHING; } marci@615: jacint@631: /// Resets the edge map of the flows to _flow. jacint@631: jacint@631: /// Resets the edge map of the flows to _flow. jacint@631: /// marci@647: void resetFlow(FlowMap& _flow) { flow=&_flow; status=AFTER_NOTHING; } marci@478: marci@478: marci@478: private: marci@478: marci@478: int push(Node w, VecStack& active) { marci@615: marci@478: int lev=level[w]; marci@478: Num exc=excess[w]; marci@478: int newlevel=n; //bound on the next level of w marci@615: marci@478: OutEdgeIt e; marci@478: for(g->first(e,w); g->valid(e); g->next(e)) { marci@615: marci@615: if ( (*flow)[e] >= (*capacity)[e] ) continue; marci@615: Node v=g->head(e); marci@615: marci@478: if( lev > level[v] ) { //Push is allowed now marci@615: marci@478: if ( excess[v]<=0 && v!=t && v!=s ) { marci@478: int lev_v=level[v]; marci@478: active[lev_v].push(v); marci@478: } marci@615: marci@478: Num cap=(*capacity)[e]; marci@478: Num flo=(*flow)[e]; marci@478: Num remcap=cap-flo; marci@615: marci@478: if ( remcap >= exc ) { //A nonsaturating push. marci@615: marci@478: flow->set(e, flo+exc); marci@478: excess.set(v, excess[v]+exc); marci@478: exc=0; marci@615: break; marci@615: marci@478: } else { //A saturating push. marci@478: flow->set(e, cap); marci@478: excess.set(v, excess[v]+remcap); marci@478: exc-=remcap; marci@478: } marci@478: } else if ( newlevel > level[v] ) newlevel = level[v]; marci@615: } //for out edges wv marci@615: marci@615: if ( exc > 0 ) { marci@478: InEdgeIt e; marci@478: for(g->first(e,w); g->valid(e); g->next(e)) { marci@615: marci@615: if( (*flow)[e] <= 0 ) continue; marci@615: Node v=g->tail(e); marci@615: marci@478: if( lev > level[v] ) { //Push is allowed now marci@615: marci@478: if ( excess[v]<=0 && v!=t && v!=s ) { marci@478: int lev_v=level[v]; marci@478: active[lev_v].push(v); marci@478: } marci@615: marci@478: Num flo=(*flow)[e]; marci@615: marci@478: if ( flo >= exc ) { //A nonsaturating push. marci@615: marci@478: flow->set(e, flo-exc); marci@478: excess.set(v, excess[v]+exc); marci@478: exc=0; marci@615: break; marci@478: } else { //A saturating push. marci@615: marci@478: excess.set(v, excess[v]+flo); marci@478: exc-=flo; marci@478: flow->set(e,0); marci@615: } marci@478: } else if ( newlevel > level[v] ) newlevel = level[v]; marci@478: } //for in edges vw marci@615: marci@478: } // if w still has excess after the out edge for cycle marci@615: marci@478: excess.set(w, exc); marci@615: marci@478: return newlevel; marci@485: } marci@478: marci@478: marci@647: void preflowPreproc(FlowEnum fe, VecStack& active, marci@615: VecNode& level_list, NNMap& left, NNMap& right) marci@602: { marci@615: std::queue bfs_queue; marci@478: marci@615: switch (fe) { jacint@631: case NO_FLOW: //flow is already set to const zero in this case marci@615: case ZERO_FLOW: marci@602: { marci@602: //Reverse_bfs from t, to find the starting level. marci@602: level.set(t,0); marci@602: bfs_queue.push(t); marci@615: marci@602: while (!bfs_queue.empty()) { marci@615: marci@615: Node v=bfs_queue.front(); marci@602: bfs_queue.pop(); marci@602: int l=level[v]+1; marci@615: marci@602: InEdgeIt e; marci@602: for(g->first(e,v); g->valid(e); g->next(e)) { marci@602: Node w=g->tail(e); marci@602: if ( level[w] == n && w != s ) { marci@602: bfs_queue.push(w); marci@602: Node first=level_list[l]; marci@602: if ( g->valid(first) ) left.set(first,w); marci@602: right.set(w,first); marci@602: level_list[l]=w; marci@602: level.set(w, l); marci@602: } marci@602: } marci@602: } marci@615: marci@602: //the starting flow marci@602: OutEdgeIt e; marci@615: for(g->first(e,s); g->valid(e); g->next(e)) marci@602: { marci@602: Num c=(*capacity)[e]; marci@602: if ( c <= 0 ) continue; marci@602: Node w=g->head(e); marci@615: if ( level[w] < n ) { marci@602: if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); marci@615: flow->set(e, c); marci@602: excess.set(w, excess[w]+c); marci@602: } marci@602: } marci@602: break; marci@602: } marci@615: marci@602: case GEN_FLOW: marci@615: case PRE_FLOW: marci@602: { marci@615: //Reverse_bfs from t in the residual graph, marci@602: //to find the starting level. marci@602: level.set(t,0); marci@602: bfs_queue.push(t); marci@615: marci@602: while (!bfs_queue.empty()) { marci@615: marci@615: Node v=bfs_queue.front(); marci@602: bfs_queue.pop(); marci@602: int l=level[v]+1; marci@615: marci@602: InEdgeIt e; marci@602: for(g->first(e,v); g->valid(e); g->next(e)) { marci@602: if ( (*capacity)[e] <= (*flow)[e] ) continue; marci@602: Node w=g->tail(e); marci@602: if ( level[w] == n && w != s ) { marci@602: bfs_queue.push(w); marci@602: Node first=level_list[l]; marci@602: if ( g->valid(first) ) left.set(first,w); marci@602: right.set(w,first); marci@602: level_list[l]=w; marci@602: level.set(w, l); marci@602: } marci@602: } marci@615: marci@602: OutEdgeIt f; marci@602: for(g->first(f,v); g->valid(f); g->next(f)) { marci@602: if ( 0 >= (*flow)[f] ) continue; marci@602: Node w=g->head(f); marci@602: if ( level[w] == n && w != s ) { marci@602: bfs_queue.push(w); marci@602: Node first=level_list[l]; marci@602: if ( g->valid(first) ) left.set(first,w); marci@602: right.set(w,first); marci@602: level_list[l]=w; marci@602: level.set(w, l); marci@602: } marci@602: } marci@602: } marci@615: marci@615: marci@602: //the starting flow marci@602: OutEdgeIt e; marci@615: for(g->first(e,s); g->valid(e); g->next(e)) marci@602: { marci@602: Num rem=(*capacity)[e]-(*flow)[e]; marci@602: if ( rem <= 0 ) continue; marci@602: Node w=g->head(e); marci@615: if ( level[w] < n ) { marci@602: if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); marci@615: flow->set(e, (*capacity)[e]); marci@602: excess.set(w, excess[w]+rem); marci@602: } marci@602: } marci@615: marci@602: InEdgeIt f; marci@615: for(g->first(f,s); g->valid(f); g->next(f)) marci@602: { marci@602: if ( (*flow)[f] <= 0 ) continue; marci@602: Node w=g->tail(f); marci@615: if ( level[w] < n ) { marci@602: if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); marci@602: excess.set(w, excess[w]+(*flow)[f]); marci@615: flow->set(f, 0); marci@602: } marci@615: } marci@602: break; marci@615: } //case PRE_FLOW marci@602: } marci@602: } //preflowPreproc marci@478: marci@478: marci@478: marci@615: void relabel(Node w, int newlevel, VecStack& active, marci@615: VecNode& level_list, NNMap& left, marci@615: NNMap& right, int& b, int& k, bool what_heur ) marci@478: { marci@478: marci@615: Num lev=level[w]; marci@615: marci@478: Node right_n=right[w]; marci@478: Node left_n=left[w]; marci@615: marci@478: //unlacing starts marci@478: if ( g->valid(right_n) ) { marci@478: if ( g->valid(left_n) ) { marci@478: right.set(left_n, right_n); marci@478: left.set(right_n, left_n); marci@478: } else { marci@615: level_list[lev]=right_n; marci@478: left.set(right_n, INVALID); marci@615: } marci@478: } else { marci@478: if ( g->valid(left_n) ) { marci@478: right.set(left_n, INVALID); marci@615: } else { marci@615: level_list[lev]=INVALID; marci@615: } marci@615: } marci@478: //unlacing ends marci@615: marci@478: if ( !g->valid(level_list[lev]) ) { marci@615: marci@478: //gapping starts marci@478: for (int i=lev; i!=k ; ) { marci@478: Node v=level_list[++i]; marci@478: while ( g->valid(v) ) { marci@478: level.set(v,n); marci@478: v=right[v]; marci@478: } marci@478: level_list[i]=INVALID; marci@478: if ( !what_heur ) { marci@478: while ( !active[i].empty() ) { marci@478: active[i].pop(); //FIXME: ezt szebben kene marci@478: } marci@615: } marci@478: } marci@615: marci@478: level.set(w,n); marci@478: b=lev-1; marci@478: k=b; marci@478: //gapping ends marci@615: marci@478: } else { marci@615: marci@615: if ( newlevel == n ) level.set(w,n); marci@478: else { marci@478: level.set(w,++newlevel); marci@478: active[newlevel].push(w); marci@478: if ( what_heur ) b=newlevel; marci@478: if ( k < newlevel ) ++k; //now k=newlevel marci@478: Node first=level_list[newlevel]; marci@478: if ( g->valid(first) ) left.set(first,w); marci@478: right.set(w,first); marci@478: left.set(w,INVALID); marci@478: level_list[newlevel]=w; marci@478: } marci@478: } marci@615: marci@478: } //relabel marci@478: marci@478: marci@615: template marci@478: class DistanceMap { marci@478: protected: marci@478: const MapGraphWrapper* g; marci@615: typename MapGraphWrapper::template NodeMap dist; marci@478: public: marci@478: DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { } marci@615: void set(const typename MapGraphWrapper::Node& n, int a) { marci@615: dist.set(n, a); marci@478: } marci@647: int operator[](const typename MapGraphWrapper::Node& n) const { marci@647: return dist[n]; marci@647: } marci@615: // int get(const typename MapGraphWrapper::Node& n) const { marci@485: // return dist[n]; } marci@615: // bool get(const typename MapGraphWrapper::Edge& e) const { marci@485: // return (dist.get(g->tail(e))head(e))); } marci@615: bool operator[](const typename MapGraphWrapper::Edge& e) const { marci@615: return (dist[g->tail(e)]head(e)]); marci@478: } marci@478: }; marci@615: marci@478: }; marci@478: marci@478: marci@478: template marci@647: void MaxFlow::preflowPhase1(FlowEnum fe) marci@478: { marci@615: marci@615: int heur0=(int)(H0*n); //time while running 'bound decrease' marci@485: int heur1=(int)(H1*n); //time while running 'highest label' marci@485: int heur=heur1; //starting time interval (#of relabels) marci@485: int numrelabel=0; marci@615: marci@615: bool what_heur=1; marci@485: //It is 0 in case 'bound decrease' and 1 in case 'highest label' marci@478: marci@615: bool end=false; marci@615: //Needed for 'bound decrease', true means no active nodes are above bound marci@615: //b. marci@478: marci@485: int k=n-2; //bound on the highest level under n containing a node marci@485: int b=k; //bound on the highest level under n of an active node marci@615: marci@485: VecStack active(n); marci@615: marci@485: NNMap left(*g, INVALID); marci@485: NNMap right(*g, INVALID); marci@485: VecNode level_list(n,INVALID); marci@485: //List of the nodes in level ifirst(v); g->valid(v); g->next(v)) level.set(v,n); marci@485: //setting each node to level n marci@615: jacint@631: if ( fe == NO_FLOW ) { jacint@631: EdgeIt e; jacint@631: for(g->first(e); g->valid(e); g->next(e)) flow->set(e,0); jacint@631: } jacint@631: jacint@631: switch (fe) { //computing the excess marci@615: case PRE_FLOW: marci@485: { marci@485: NodeIt v; marci@485: for(g->first(v); g->valid(v); g->next(v)) { marci@478: Num exc=0; marci@615: marci@478: InEdgeIt e; marci@485: for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e]; marci@478: OutEdgeIt f; marci@485: for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f]; marci@615: marci@615: excess.set(v,exc); marci@615: marci@485: //putting the active nodes into the stack marci@485: int lev=level[v]; marci@485: if ( exc > 0 && lev < n && v != t ) active[lev].push(v); marci@478: } marci@478: break; marci@478: } marci@485: case GEN_FLOW: marci@485: { jacint@631: NodeIt v; jacint@631: for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0); jacint@631: marci@485: Num exc=0; marci@485: InEdgeIt e; marci@485: for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e]; marci@485: OutEdgeIt f; marci@485: for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f]; marci@615: excess.set(t,exc); marci@485: break; marci@485: } jacint@631: case ZERO_FLOW: jacint@631: case NO_FLOW: jacint@631: { jacint@631: NodeIt v; jacint@631: for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0); jacint@631: break; jacint@631: } marci@485: } marci@615: marci@615: preflowPreproc(fe, active, level_list, left, right); marci@615: //End of preprocessing marci@615: marci@615: marci@485: //Push/relabel on the highest level active nodes. marci@485: while ( true ) { marci@485: if ( b == 0 ) { marci@485: if ( !what_heur && !end && k > 0 ) { marci@485: b=k; marci@485: end=true; marci@485: } else break; marci@485: } marci@615: marci@615: if ( active[b].empty() ) --b; marci@485: else { marci@615: end=false; marci@485: Node w=active[b].top(); marci@485: active[b].pop(); marci@485: int newlevel=push(w,active); marci@615: if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list, marci@485: left, right, b, k, what_heur); marci@615: marci@615: ++numrelabel; marci@485: if ( numrelabel >= heur ) { marci@485: numrelabel=0; marci@485: if ( what_heur ) { marci@485: what_heur=0; marci@485: heur=heur0; marci@485: end=false; marci@485: } else { marci@485: what_heur=1; marci@485: heur=heur1; marci@615: b=k; marci@485: } marci@478: } marci@615: } marci@615: } marci@647: marci@647: status=AFTER_PRE_FLOW_PHASE_1; marci@485: } marci@478: marci@478: marci@478: marci@478: template jacint@631: void MaxFlow::preflowPhase2() marci@478: { marci@615: marci@485: int k=n-2; //bound on the highest level under n containing a node marci@485: int b=k; //bound on the highest level under n of an active node marci@615: marci@485: VecStack active(n); marci@485: level.set(s,0); marci@485: std::queue bfs_queue; marci@485: bfs_queue.push(s); marci@615: marci@485: while (!bfs_queue.empty()) { marci@615: marci@615: Node v=bfs_queue.front(); marci@485: bfs_queue.pop(); marci@485: int l=level[v]+1; marci@615: marci@485: InEdgeIt e; marci@485: for(g->first(e,v); g->valid(e); g->next(e)) { marci@485: if ( (*capacity)[e] <= (*flow)[e] ) continue; marci@485: Node u=g->tail(e); marci@615: if ( level[u] >= n ) { marci@485: bfs_queue.push(u); marci@485: level.set(u, l); marci@485: if ( excess[u] > 0 ) active[l].push(u); marci@478: } marci@478: } marci@615: marci@485: OutEdgeIt f; marci@485: for(g->first(f,v); g->valid(f); g->next(f)) { marci@485: if ( 0 >= (*flow)[f] ) continue; marci@485: Node u=g->head(f); marci@615: if ( level[u] >= n ) { marci@485: bfs_queue.push(u); marci@485: level.set(u, l); marci@485: if ( excess[u] > 0 ) active[l].push(u); marci@485: } marci@485: } marci@485: } marci@485: b=n-2; marci@478: marci@485: while ( true ) { marci@615: marci@485: if ( b == 0 ) break; marci@478: marci@615: if ( active[b].empty() ) --b; marci@485: else { marci@485: Node w=active[b].top(); marci@485: active[b].pop(); marci@615: int newlevel=push(w,active); marci@478: marci@485: //relabel marci@485: if ( excess[w] > 0 ) { marci@485: level.set(w,++newlevel); marci@485: active[newlevel].push(w); marci@485: b=newlevel; marci@485: } marci@485: } // if stack[b] is nonempty marci@485: } // while(true) marci@647: marci@647: status=AFTER_PRE_FLOW_PHASE_2; marci@485: } marci@478: marci@478: marci@478: marci@478: template marci@615: bool MaxFlow::augmentOnShortestPath() marci@478: { marci@485: ResGW res_graph(*g, *capacity, *flow); marci@485: bool _augment=false; marci@615: marci@485: //ReachedMap level(res_graph); marci@485: FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); marci@485: BfsIterator bfs(res_graph, level); marci@485: bfs.pushAndSetReached(s); marci@615: marci@615: typename ResGW::template NodeMap pred(res_graph); marci@485: pred.set(s, INVALID); marci@615: marci@485: typename ResGW::template NodeMap free(res_graph); marci@615: marci@485: //searching for augmenting path marci@615: while ( !bfs.finished() ) { marci@485: ResGWOutEdgeIt e=bfs; marci@485: if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) { marci@485: Node v=res_graph.tail(e); marci@485: Node w=res_graph.head(e); marci@485: pred.set(w, e); marci@485: if (res_graph.valid(pred[v])) { marci@485: free.set(w, std::min(free[v], res_graph.resCap(e))); marci@485: } else { marci@615: free.set(w, res_graph.resCap(e)); marci@478: } marci@485: if (res_graph.head(e)==t) { _augment=true; break; } marci@485: } marci@615: marci@485: ++bfs; marci@485: } //end of searching augmenting path marci@478: marci@485: if (_augment) { marci@485: Node n=t; marci@485: Num augment_value=free[t]; marci@615: while (res_graph.valid(pred[n])) { marci@485: ResGWEdge e=pred[n]; marci@615: res_graph.augment(e, augment_value); marci@485: n=res_graph.tail(e); marci@478: } marci@485: } marci@478: marci@647: status=AFTER_AUGMENTING; marci@485: return _augment; marci@485: } marci@478: marci@478: marci@647: template marci@647: bool MaxFlow::augmentOnShortestPath2() marci@647: { marci@647: ResGW res_graph(*g, *capacity, *flow); marci@647: bool _augment=false; marci@478: marci@656: if (status!=AFTER_FAST_AUGMENTING) { marci@656: FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); marci@656: number_of_augmentations=1; marci@647: } else { marci@647: ++number_of_augmentations; marci@647: } marci@647: TrickyReachedMap marci@647: tricky_reached_map(level, number_of_augmentations); marci@647: //ReachedMap level(res_graph); marci@647: // FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); marci@647: BfsIterator > marci@647: bfs(res_graph, tricky_reached_map); marci@647: bfs.pushAndSetReached(s); marci@478: marci@647: typename ResGW::template NodeMap pred(res_graph); marci@647: pred.set(s, INVALID); marci@478: marci@647: typename ResGW::template NodeMap free(res_graph); marci@647: marci@647: //searching for augmenting path marci@647: while ( !bfs.finished() ) { marci@647: ResGWOutEdgeIt e=bfs; marci@647: if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) { marci@647: Node v=res_graph.tail(e); marci@647: Node w=res_graph.head(e); marci@647: pred.set(w, e); marci@647: if (res_graph.valid(pred[v])) { marci@647: free.set(w, std::min(free[v], res_graph.resCap(e))); marci@647: } else { marci@647: free.set(w, res_graph.resCap(e)); marci@647: } marci@647: if (res_graph.head(e)==t) { _augment=true; break; } marci@647: } marci@647: marci@647: ++bfs; marci@647: } //end of searching augmenting path marci@647: marci@647: if (_augment) { marci@647: Node n=t; marci@647: Num augment_value=free[t]; marci@647: while (res_graph.valid(pred[n])) { marci@647: ResGWEdge e=pred[n]; marci@647: res_graph.augment(e, augment_value); marci@647: n=res_graph.tail(e); marci@647: } marci@647: } marci@647: marci@656: status=AFTER_FAST_AUGMENTING; marci@647: return _augment; marci@647: } marci@478: marci@478: marci@478: template marci@615: template marci@615: bool MaxFlow::augmentOnBlockingFlow() marci@615: { marci@485: typedef MutableGraph MG; marci@485: bool _augment=false; marci@478: marci@485: ResGW res_graph(*g, *capacity, *flow); marci@478: marci@485: //bfs for distances on the residual graph marci@485: //ReachedMap level(res_graph); marci@485: FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); marci@485: BfsIterator bfs(res_graph, level); marci@485: bfs.pushAndSetReached(s); marci@615: typename ResGW::template NodeMap marci@485: dist(res_graph); //filled up with 0's marci@478: marci@485: //F will contain the physical copy of the residual graph marci@485: //with the set of edges which are on shortest paths marci@485: MG F; marci@615: typename ResGW::template NodeMap marci@485: res_graph_to_F(res_graph); marci@485: { marci@485: typename ResGW::NodeIt n; marci@485: for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) { marci@485: res_graph_to_F.set(n, F.addNode()); marci@478: } marci@485: } marci@478: marci@485: typename MG::Node sF=res_graph_to_F[s]; marci@485: typename MG::Node tF=res_graph_to_F[t]; marci@485: typename MG::template EdgeMap original_edge(F); marci@485: typename MG::template EdgeMap residual_capacity(F); marci@478: marci@615: while ( !bfs.finished() ) { marci@485: ResGWOutEdgeIt e=bfs; marci@485: if (res_graph.valid(e)) { marci@485: if (bfs.isBNodeNewlyReached()) { marci@485: dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1); marci@615: typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], marci@615: res_graph_to_F[res_graph.head(e)]); marci@485: original_edge.update(); marci@485: original_edge.set(f, e); marci@485: residual_capacity.update(); marci@485: residual_capacity.set(f, res_graph.resCap(e)); marci@485: } else { marci@485: if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) { marci@615: typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], marci@615: res_graph_to_F[res_graph.head(e)]); marci@478: original_edge.update(); marci@478: original_edge.set(f, e); marci@478: residual_capacity.update(); marci@478: residual_capacity.set(f, res_graph.resCap(e)); marci@478: } marci@478: } marci@485: } marci@485: ++bfs; marci@485: } //computing distances from s in the residual graph marci@478: marci@485: bool __augment=true; marci@478: marci@485: while (__augment) { marci@485: __augment=false; marci@485: //computing blocking flow with dfs marci@485: DfsIterator< MG, typename MG::template NodeMap > dfs(F); marci@485: typename MG::template NodeMap pred(F); marci@485: pred.set(sF, INVALID); marci@485: //invalid iterators for sources marci@478: marci@485: typename MG::template NodeMap free(F); marci@478: marci@615: dfs.pushAndSetReached(sF); marci@485: while (!dfs.finished()) { marci@485: ++dfs; marci@485: if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) { marci@485: if (dfs.isBNodeNewlyReached()) { marci@485: typename MG::Node v=F.aNode(dfs); marci@485: typename MG::Node w=F.bNode(dfs); marci@485: pred.set(w, dfs); marci@485: if (F.valid(pred[v])) { marci@485: free.set(w, std::min(free[v], residual_capacity[dfs])); marci@485: } else { marci@615: free.set(w, residual_capacity[dfs]); marci@485: } marci@615: if (w==tF) { marci@615: __augment=true; marci@485: _augment=true; marci@615: break; marci@485: } marci@615: marci@485: } else { marci@485: F.erase(/*typename MG::OutEdgeIt*/(dfs)); marci@485: } marci@615: } marci@485: } marci@485: marci@485: if (__augment) { marci@485: typename MG::Node n=tF; marci@485: Num augment_value=free[tF]; marci@615: while (F.valid(pred[n])) { marci@485: typename MG::Edge e=pred[n]; marci@615: res_graph.augment(original_edge[e], augment_value); marci@485: n=F.tail(e); marci@615: if (residual_capacity[e]==augment_value) marci@615: F.erase(e); marci@615: else marci@485: residual_capacity.set(e, residual_capacity[e]-augment_value); marci@478: } marci@485: } marci@615: marci@485: } marci@615: marci@647: status=AFTER_AUGMENTING; marci@485: return _augment; marci@485: } marci@478: marci@478: marci@478: marci@478: marci@478: template marci@615: bool MaxFlow::augmentOnBlockingFlow2() marci@478: { marci@485: bool _augment=false; marci@478: marci@485: ResGW res_graph(*g, *capacity, *flow); marci@615: marci@485: //ReachedMap level(res_graph); marci@485: FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); marci@485: BfsIterator bfs(res_graph, level); marci@478: marci@485: bfs.pushAndSetReached(s); marci@485: DistanceMap dist(res_graph); marci@615: while ( !bfs.finished() ) { marci@485: ResGWOutEdgeIt e=bfs; marci@485: if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) { marci@485: dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1); marci@485: } marci@485: ++bfs; marci@485: } //computing distances from s in the residual graph marci@478: marci@478: //Subgraph containing the edges on some shortest paths marci@485: ConstMap true_map(true); marci@615: typedef SubGraphWrapper, marci@485: DistanceMap > FilterResGW; marci@485: FilterResGW filter_res_graph(res_graph, true_map, dist); marci@478: marci@615: //Subgraph, which is able to delete edges which are already marci@485: //met by the dfs marci@615: typename FilterResGW::template NodeMap marci@485: first_out_edges(filter_res_graph); marci@485: typename FilterResGW::NodeIt v; marci@615: for(filter_res_graph.first(v); filter_res_graph.valid(v); marci@615: filter_res_graph.next(v)) marci@478: { marci@478: typename FilterResGW::OutEdgeIt e; marci@478: filter_res_graph.first(e, v); marci@478: first_out_edges.set(v, e); marci@478: } marci@485: typedef ErasingFirstGraphWrapper > ErasingResGW; marci@485: ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges); marci@478: marci@485: bool __augment=true; marci@478: marci@485: while (__augment) { marci@478: marci@485: __augment=false; marci@485: //computing blocking flow with dfs marci@615: DfsIterator< ErasingResGW, marci@615: typename ErasingResGW::template NodeMap > marci@485: dfs(erasing_res_graph); marci@485: typename ErasingResGW:: marci@615: template NodeMap marci@615: pred(erasing_res_graph); marci@485: pred.set(s, INVALID); marci@485: //invalid iterators for sources marci@478: marci@615: typename ErasingResGW::template NodeMap marci@485: free1(erasing_res_graph); marci@478: marci@615: dfs.pushAndSetReached marci@615: ///\bug hugo 0.2 marci@615: (typename ErasingResGW::Node marci@615: (typename FilterResGW::Node marci@615: (typename ResGW::Node(s) marci@615: ) marci@615: ) marci@615: ); marci@485: while (!dfs.finished()) { marci@485: ++dfs; marci@615: if (erasing_res_graph.valid(typename ErasingResGW::OutEdgeIt(dfs))) marci@615: { marci@478: if (dfs.isBNodeNewlyReached()) { marci@615: marci@478: typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs); marci@478: typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs); marci@478: marci@478: pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs)); marci@478: if (erasing_res_graph.valid(pred[v])) { marci@615: free1.set marci@615: (w, std::min(free1[v], res_graph.resCap marci@615: (typename ErasingResGW::OutEdgeIt(dfs)))); marci@478: } else { marci@615: free1.set marci@615: (w, res_graph.resCap marci@615: (typename ErasingResGW::OutEdgeIt(dfs))); marci@478: } marci@615: marci@615: if (w==t) { marci@615: __augment=true; marci@478: _augment=true; marci@615: break; marci@478: } marci@478: } else { marci@478: erasing_res_graph.erase(dfs); marci@478: } marci@478: } marci@615: } marci@478: marci@485: if (__augment) { marci@615: typename ErasingResGW::Node marci@615: n=typename FilterResGW::Node(typename ResGW::Node(t)); marci@485: // typename ResGW::NodeMap a(res_graph); marci@485: // typename ResGW::Node b; marci@485: // Num j=a[b]; marci@485: // typename FilterResGW::NodeMap a1(filter_res_graph); marci@485: // typename FilterResGW::Node b1; marci@485: // Num j1=a1[b1]; marci@485: // typename ErasingResGW::NodeMap a2(erasing_res_graph); marci@485: // typename ErasingResGW::Node b2; marci@485: // Num j2=a2[b2]; marci@485: Num augment_value=free1[n]; marci@615: while (erasing_res_graph.valid(pred[n])) { marci@485: typename ErasingResGW::OutEdgeIt e=pred[n]; marci@485: res_graph.augment(e, augment_value); marci@485: n=erasing_res_graph.tail(e); marci@485: if (res_graph.resCap(e)==0) marci@485: erasing_res_graph.erase(e); marci@478: } marci@478: } marci@615: marci@615: } //while (__augment) marci@615: marci@647: status=AFTER_AUGMENTING; marci@485: return _augment; marci@485: } marci@478: marci@478: marci@478: } //namespace hugo marci@478: marci@480: #endif //HUGO_MAX_FLOW_H marci@478: marci@478: marci@478: marci@478: