jacint@102: // -*- C++ -*- jacint@102: /* jacint@109: preflow_h5.h jacint@102: by jacint. jacint@109: Heuristics: jacint@109: 2 phase jacint@109: gap jacint@109: list 'level_list' on the nodes on level i implemented by hand jacint@109: highest label jacint@109: relevel: in phase 0, after BFS*n relabels, it runs a reverse jacint@109: bfs from t in the res graph to relevel the nodes reachable from t. jacint@109: BFS is initialized to 20 jacint@102: jacint@109: Due to the last heuristic, this algorithm is quite fast on very jacint@109: sparse graphs, but relatively bad on even the dense graphs. jacint@102: jacint@109: 'NodeMap cut' is a member, in this way we can count it fast, after jacint@109: the algorithm was run. jacint@102: jacint@109: The constructor runs the algorithm. jacint@102: jacint@109: Members: jacint@102: jacint@109: T maxFlow() : returns the value of a maximum flow jacint@102: jacint@109: T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e) jacint@102: jacint@109: FlowMap Flow() : returns the fixed maximum flow x jacint@102: jacint@109: void Flow(FlowMap& _flow ) : returns the fixed maximum flow x jacint@102: jacint@109: void minMinCut(CutMap& M) : sets M to the characteristic vector of the jacint@102: minimum min cut. M should be a map of bools initialized to false. jacint@102: jacint@109: void maxMinCut(CutMap& M) : sets M to the characteristic vector of the jacint@102: maximum min cut. M should be a map of bools initialized to false. jacint@102: jacint@109: void minCut(CutMap& M) : fast function, sets M to the characteristic jacint@109: vector of a minimum cut. jacint@109: jacint@109: Different member from the other preflow_hl-s (here we have a member jacint@109: 'NodeMap cut'). jacint@109: jacint@109: CutMap minCut() : fast function, giving the characteristic jacint@109: vector of a minimum cut. jacint@109: jacint@109: jacint@102: */ jacint@102: jacint@102: #ifndef PREFLOW_HL4_H jacint@102: #define PREFLOW_HL4_H jacint@102: jacint@109: #define BFS 20 jacint@109: jacint@102: #include jacint@102: #include jacint@102: jacint@109: #include //for test jacint@109: alpar@105: namespace hugo { jacint@102: jacint@102: template , jacint@109: typename CutMap=typename Graph::NodeMap, jacint@102: typename CapMap=typename Graph::EdgeMap > jacint@102: class preflow_hl4 { jacint@102: jacint@102: typedef typename Graph::NodeIt NodeIt; jacint@102: typedef typename Graph::EdgeIt EdgeIt; jacint@102: typedef typename Graph::EachNodeIt EachNodeIt; jacint@102: typedef typename Graph::OutEdgeIt OutEdgeIt; jacint@102: typedef typename Graph::InEdgeIt InEdgeIt; jacint@102: jacint@102: Graph& G; jacint@102: NodeIt s; jacint@102: NodeIt t; jacint@102: FlowMap flow; jacint@102: CapMap& capacity; jacint@109: CutMap cut; jacint@102: T value; jacint@102: jacint@102: public: jacint@102: jacint@109: double time; jacint@109: jacint@102: preflow_hl4(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) : jacint@109: G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity), jacint@109: cut(G, false) { jacint@102: jacint@109: bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd jacint@102: int n=G.nodeNum(); jacint@109: int relabel=0; jacint@109: int heur=(int)BFS*n; jacint@102: int k=n-2; jacint@102: int b=k; jacint@102: /* jacint@102: b is a bound on the highest level of the stack. jacint@102: k is a bound on the highest nonempty level i < n. jacint@102: */ jacint@102: jacint@102: typename Graph::NodeMap level(G,n); jacint@102: typename Graph::NodeMap excess(G); jacint@109: jacint@109: std::vector active(n); jacint@109: typename Graph::NodeMap next(G); jacint@102: //Stack of the active nodes in level i < n. jacint@102: //We use it in both phases. jacint@102: jacint@102: typename Graph::NodeMap left(G); jacint@102: typename Graph::NodeMap right(G); jacint@102: std::vector level_list(n); jacint@102: /* jacint@102: Needed for the list of the nodes in level i. jacint@102: */ jacint@102: jacint@102: /*Reverse_bfs from t, to find the starting level.*/ jacint@102: level.set(t,0); jacint@102: std::queue bfs_queue; jacint@102: bfs_queue.push(t); jacint@102: jacint@102: while (!bfs_queue.empty()) { jacint@102: jacint@102: NodeIt v=bfs_queue.front(); jacint@102: bfs_queue.pop(); jacint@102: int l=level.get(v)+1; jacint@102: jacint@102: for(InEdgeIt e=G.template first(v); e.valid(); ++e) { jacint@102: NodeIt w=G.tail(e); jacint@109: if ( level.get(w) == n && w !=s ) { jacint@102: bfs_queue.push(w); jacint@102: NodeIt first=level_list[l]; jacint@102: if ( first != 0 ) left.set(first,w); jacint@102: right.set(w,first); jacint@102: level_list[l]=w; jacint@102: level.set(w, l); jacint@102: } jacint@102: } jacint@102: } jacint@102: jacint@102: level.set(s,n); jacint@102: jacint@102: jacint@102: /* Starting flow. It is everywhere 0 at the moment. */ jacint@102: for(OutEdgeIt e=G.template first(s); e.valid(); ++e) jacint@102: { jacint@102: T c=capacity.get(e); jacint@102: if ( c == 0 ) continue; jacint@102: NodeIt w=G.head(e); jacint@102: if ( level.get(w) < n ) { jacint@109: if ( excess.get(w) == 0 && w!=t ) { jacint@109: next.set(w,active[level.get(w)]); jacint@109: active[level.get(w)]=w; jacint@109: } jacint@102: flow.set(e, c); jacint@102: excess.set(w, excess.get(w)+c); jacint@102: } jacint@102: } jacint@102: /* jacint@102: End of preprocessing jacint@102: */ jacint@102: jacint@102: jacint@102: /* jacint@102: Push/relabel on the highest level active nodes. jacint@102: */ jacint@102: while ( true ) { jacint@102: jacint@102: if ( b == 0 ) { jacint@102: if ( phase ) break; jacint@102: jacint@102: /* jacint@102: In the end of phase 0 we apply a bfs from s in jacint@102: the residual graph. jacint@102: */ jacint@102: phase=1; jacint@109: jacint@109: //Now have a min cut. jacint@109: for( EachNodeIt v=G.template first(); jacint@109: v.valid(); ++v) jacint@109: if (level.get(v) >= n ) cut.set(v,true); jacint@109: jacint@109: time=currTime(); jacint@102: level.set(s,0); jacint@102: std::queue bfs_queue; jacint@102: bfs_queue.push(s); jacint@102: jacint@102: while (!bfs_queue.empty()) { jacint@102: jacint@102: NodeIt v=bfs_queue.front(); jacint@102: bfs_queue.pop(); jacint@102: int l=level.get(v)+1; jacint@102: jacint@102: for(InEdgeIt e=G.template first(v); e.valid(); ++e) { jacint@102: if ( capacity.get(e) == flow.get(e) ) continue; jacint@102: NodeIt u=G.tail(e); jacint@102: if ( level.get(u) >= n ) { jacint@102: bfs_queue.push(u); jacint@102: level.set(u, l); jacint@109: if ( excess.get(u) > 0 ) { jacint@109: next.set(u,active[l]); jacint@109: active[l]=u; jacint@109: } jacint@102: } jacint@102: } jacint@102: jacint@102: for(OutEdgeIt e=G.template first(v); e.valid(); ++e) { jacint@102: if ( 0 == flow.get(e) ) continue; jacint@102: NodeIt u=G.head(e); jacint@102: if ( level.get(u) >= n ) { jacint@102: bfs_queue.push(u); jacint@102: level.set(u, l); jacint@109: if ( excess.get(u) > 0 ) { jacint@109: next.set(u,active[l]); jacint@109: active[l]=u; jacint@109: } jacint@102: } jacint@102: } jacint@102: } jacint@102: b=n-2; jacint@102: } jacint@102: jacint@102: jacint@109: if ( active[b] == 0 ) --b; jacint@102: else { jacint@102: jacint@109: NodeIt w=active[b]; jacint@109: active[b]=next.get(w); jacint@102: int lev=level.get(w); jacint@102: T exc=excess.get(w); jacint@109: int newlevel=n; //bound on the next level of w. jacint@102: jacint@102: for(OutEdgeIt e=G.template first(w); e.valid(); ++e) { jacint@102: jacint@102: if ( flow.get(e) == capacity.get(e) ) continue; jacint@102: NodeIt v=G.head(e); jacint@102: //e=wv jacint@102: jacint@102: if( lev > level.get(v) ) { jacint@102: /*Push is allowed now*/ jacint@102: jacint@109: if ( excess.get(v)==0 && v!=t && v!=s ) { jacint@109: int lev_v=level.get(v); jacint@109: next.set(v,active[lev_v]); jacint@109: active[lev_v]=v; jacint@109: } jacint@102: jacint@102: T cap=capacity.get(e); jacint@102: T flo=flow.get(e); jacint@102: T remcap=cap-flo; jacint@102: jacint@102: if ( remcap >= exc ) { jacint@102: /*A nonsaturating push.*/ jacint@102: jacint@102: flow.set(e, flo+exc); jacint@102: excess.set(v, excess.get(v)+exc); jacint@102: exc=0; jacint@102: break; jacint@102: jacint@102: } else { jacint@102: /*A saturating push.*/ jacint@102: jacint@102: flow.set(e, cap); jacint@102: excess.set(v, excess.get(v)+remcap); jacint@102: exc-=remcap; jacint@102: } jacint@102: } else if ( newlevel > level.get(v) ){ jacint@102: newlevel = level.get(v); jacint@102: } jacint@102: jacint@102: } //for out edges wv jacint@102: jacint@102: jacint@102: if ( exc > 0 ) { jacint@102: for( InEdgeIt e=G.template first(w); e.valid(); ++e) { jacint@102: jacint@102: if( flow.get(e) == 0 ) continue; jacint@102: NodeIt v=G.tail(e); jacint@102: //e=vw jacint@102: jacint@102: if( lev > level.get(v) ) { jacint@102: /*Push is allowed now*/ jacint@102: jacint@109: if ( excess.get(v)==0 && v!=t && v!=s ) { jacint@109: int lev_v=level.get(v); jacint@109: next.set(v,active[lev_v]); jacint@109: active[lev_v]=v; jacint@109: } jacint@102: jacint@102: T flo=flow.get(e); jacint@102: jacint@102: if ( flo >= exc ) { jacint@102: /*A nonsaturating push.*/ jacint@102: jacint@102: flow.set(e, flo-exc); jacint@102: excess.set(v, excess.get(v)+exc); jacint@102: exc=0; jacint@102: break; jacint@102: } else { jacint@102: /*A saturating push.*/ jacint@102: jacint@102: excess.set(v, excess.get(v)+flo); jacint@102: exc-=flo; jacint@102: flow.set(e,0); jacint@102: } jacint@102: } else if ( newlevel > level.get(v) ) { jacint@102: newlevel = level.get(v); jacint@102: } jacint@102: } //for in edges vw jacint@102: jacint@102: } // if w still has excess after the out edge for cycle jacint@102: jacint@102: excess.set(w, exc); jacint@102: jacint@102: /* jacint@102: Relabel jacint@102: */ jacint@102: jacint@102: if ( exc > 0 ) { jacint@102: //now 'lev' is the old level of w jacint@102: jacint@102: if ( phase ) { jacint@102: level.set(w,++newlevel); jacint@109: next.set(w,active[newlevel]); jacint@109: active[newlevel]=w; jacint@102: b=newlevel; jacint@102: } else { jacint@102: //unlacing jacint@102: NodeIt right_n=right.get(w); jacint@102: NodeIt left_n=left.get(w); jacint@102: jacint@102: if ( right_n != 0 ) { jacint@102: if ( left_n != 0 ) { jacint@102: right.set(left_n, right_n); jacint@102: left.set(right_n, left_n); jacint@102: } else { jacint@102: level_list[lev]=right_n; jacint@102: left.set(right_n, 0); jacint@102: } jacint@102: } else { jacint@102: if ( left_n != 0 ) { jacint@102: right.set(left_n, 0); jacint@102: } else { jacint@102: level_list[lev]=0; jacint@102: } jacint@102: } jacint@109: //unlacing ends jacint@102: jacint@102: jacint@102: if ( level_list[lev]==0 ) { jacint@102: jacint@102: for (int i=lev; i!=k ; ) { jacint@102: NodeIt v=level_list[++i]; jacint@102: while ( v != 0 ) { jacint@102: level.set(v,n); jacint@102: v=right.get(v); jacint@102: } jacint@102: level_list[i]=0; jacint@102: } jacint@102: jacint@102: level.set(w,n); jacint@102: jacint@102: b=--lev; jacint@102: k=b; jacint@102: jacint@102: } else { jacint@102: jacint@102: if ( newlevel == n ) { jacint@102: level.set(w,n); jacint@102: } else { jacint@102: jacint@102: level.set(w,++newlevel); jacint@109: next.set(w,active[newlevel]); jacint@109: active[newlevel]=w; jacint@102: b=newlevel; jacint@102: if ( k < newlevel ) ++k; jacint@102: NodeIt first=level_list[newlevel]; jacint@102: if ( first != 0 ) left.set(first,w); jacint@102: right.set(w,first); jacint@102: left.set(w,0); jacint@102: level_list[newlevel]=w; jacint@102: } jacint@102: } jacint@109: jacint@109: ++relabel; jacint@109: if ( relabel >= heur ) { jacint@109: relabel=0; jacint@109: b=n-2; jacint@109: k=b; jacint@109: jacint@109: for ( int i=1; i!=n; ++i ) { jacint@109: active[i]=0; jacint@109: level_list[i]=0; jacint@109: } jacint@109: jacint@109: //bfs from t in the res graph to relevel the nodes jacint@109: for( EachNodeIt v=G.template first(); jacint@109: v.valid(); ++v) level.set(v,n); jacint@109: jacint@109: level.set(t,0); jacint@109: std::queue bfs_queue; jacint@109: bfs_queue.push(t); jacint@109: jacint@109: while (!bfs_queue.empty()) { jacint@109: jacint@109: NodeIt v=bfs_queue.front(); jacint@109: bfs_queue.pop(); jacint@109: int l=level.get(v)+1; jacint@109: jacint@109: for(InEdgeIt e=G.template first(v); jacint@109: e.valid(); ++e) { jacint@109: if ( capacity.get(e) == flow.get(e) ) continue; jacint@109: NodeIt u=G.tail(e); jacint@109: if ( level.get(u) == n ) { jacint@109: bfs_queue.push(u); jacint@109: level.set(u, l); jacint@109: if ( excess.get(u) > 0 ) { jacint@109: next.set(u,active[l]); jacint@109: active[l]=u; jacint@109: } jacint@109: NodeIt first=level_list[l]; jacint@109: if ( first != 0 ) left.set(first,w); jacint@109: right.set(w,first); jacint@109: left.set(w,0); jacint@109: level_list[l]=w; jacint@109: } jacint@109: } jacint@109: jacint@109: jacint@109: for(OutEdgeIt e=G.template first(v); jacint@109: e.valid(); ++e) { jacint@109: if ( 0 == flow.get(e) ) continue; jacint@109: NodeIt u=G.head(e); jacint@109: if ( level.get(u) == n ) { jacint@109: bfs_queue.push(u); jacint@109: level.set(u, l); jacint@109: if ( excess.get(u) > 0 ) { jacint@109: next.set(u,active[l]); jacint@109: active[l]=u; jacint@109: } jacint@109: NodeIt first=level_list[l]; jacint@109: if ( first != 0 ) left.set(first,w); jacint@109: right.set(w,first); jacint@109: left.set(w,0); jacint@109: level_list[l]=w; jacint@109: } jacint@109: } jacint@109: } jacint@109: jacint@109: level.set(s,n); jacint@109: } jacint@109: jacint@102: } //phase 0 jacint@102: } // if ( exc > 0 ) jacint@109: jacint@102: jacint@102: } // if stack[b] is nonempty jacint@102: jacint@102: } // while(true) jacint@102: jacint@102: jacint@102: value = excess.get(t); jacint@102: /*Max flow value.*/ jacint@102: jacint@102: jacint@102: } //void run() jacint@102: jacint@102: jacint@102: jacint@102: jacint@102: jacint@102: /* jacint@102: Returns the maximum value of a flow. jacint@102: */ jacint@102: jacint@109: T maxFlow() { jacint@102: return value; jacint@102: } jacint@102: jacint@102: jacint@102: jacint@102: /* jacint@102: For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e). jacint@102: */ jacint@102: jacint@109: T flowOnEdge(EdgeIt e) { jacint@102: return flow.get(e); jacint@102: } jacint@102: jacint@102: jacint@102: jacint@109: FlowMap Flow() { jacint@102: return flow; jacint@102: } jacint@102: jacint@102: jacint@102: jacint@109: void Flow(FlowMap& _flow ) { jacint@102: for(EachNodeIt v=G.template first() ; v.valid(); ++v) jacint@102: _flow.set(v,flow.get(v)); jacint@102: } jacint@102: jacint@102: jacint@102: jacint@102: /* jacint@102: Returns the minimum min cut, by a bfs from s in the residual graph. jacint@102: */ jacint@102: jacint@109: template jacint@109: void minMinCut(_CutMap& M) { jacint@102: jacint@102: std::queue queue; jacint@102: jacint@102: M.set(s,true); jacint@102: queue.push(s); jacint@102: jacint@102: while (!queue.empty()) { jacint@102: NodeIt w=queue.front(); jacint@102: queue.pop(); jacint@102: jacint@102: for(OutEdgeIt e=G.template first(w) ; e.valid(); ++e) { jacint@102: NodeIt v=G.head(e); jacint@102: if (!M.get(v) && flow.get(e) < capacity.get(e) ) { jacint@102: queue.push(v); jacint@102: M.set(v, true); jacint@102: } jacint@102: } jacint@102: jacint@102: for(InEdgeIt e=G.template first(w) ; e.valid(); ++e) { jacint@102: NodeIt v=G.tail(e); jacint@102: if (!M.get(v) && flow.get(e) > 0 ) { jacint@102: queue.push(v); jacint@102: M.set(v, true); jacint@102: } jacint@102: } jacint@102: jacint@102: } jacint@102: jacint@102: } jacint@102: jacint@102: jacint@102: jacint@102: /* jacint@102: Returns the maximum min cut, by a reverse bfs jacint@102: from t in the residual graph. jacint@102: */ jacint@102: jacint@109: template jacint@109: void maxMinCut(_CutMap& M) { jacint@102: jacint@102: std::queue queue; jacint@102: jacint@102: M.set(t,true); jacint@102: queue.push(t); jacint@102: jacint@102: while (!queue.empty()) { jacint@102: NodeIt w=queue.front(); jacint@102: queue.pop(); jacint@102: jacint@102: for(InEdgeIt e=G.template first(w) ; e.valid(); ++e) { jacint@102: NodeIt v=G.tail(e); jacint@102: if (!M.get(v) && flow.get(e) < capacity.get(e) ) { jacint@102: queue.push(v); jacint@102: M.set(v, true); jacint@102: } jacint@102: } jacint@102: jacint@102: for(OutEdgeIt e=G.template first(w) ; e.valid(); ++e) { jacint@102: NodeIt v=G.head(e); jacint@102: if (!M.get(v) && flow.get(e) > 0 ) { jacint@102: queue.push(v); jacint@102: M.set(v, true); jacint@102: } jacint@102: } jacint@102: } jacint@102: jacint@102: for(EachNodeIt v=G.template first() ; v.valid(); ++v) { jacint@102: M.set(v, !M.get(v)); jacint@102: } jacint@102: jacint@102: } jacint@102: jacint@102: jacint@109: template jacint@109: void minCut(_CutMap& M) { jacint@109: for( EachNodeIt v=G.template first(); jacint@109: v.valid(); ++v) jacint@109: M.set(v, cut.get(v)); jacint@102: } jacint@102: jacint@109: jacint@109: CutMap minCut() { jacint@109: return cut; jacint@109: } jacint@102: jacint@102: jacint@102: }; jacint@109: }//namespace marci jacint@102: #endif jacint@102: jacint@102: jacint@102: jacint@102: