jacint@140: // -*- C++ -*- jacint@140: jacint@140: //hogy kell azt megcsinalni hogy a konstruktorban a FlowMap-nek legyen default erteke (a 0 map) jacint@140: jacint@140: /* jacint@140: preflow_aug.h jacint@140: by jacint. jacint@140: jacint@140: The same as preflow.h, the only difference is that here jacint@140: we start from a given flow. jacint@140: jacint@140: jacint@140: The constructor runs the algorithm: jacint@140: jacint@140: template , jacint@140: typename Graph_EdgeMap_T_2=typename Graph::EdgeMap > jacint@140: PreflowAug(Graph G, G_NodeIt s, G_NodeIt t, G_EdgeMap_T_1 capacity, jacint@140: G_EdgeMap_T_2 flow, bool flow_type) jacint@140: jacint@140: 'capacity' must be non-negative, if flow_type=0 then jacint@140: 'flow' must be a flow, otherwise it must be a preflow. jacint@140: jacint@140: jacint@140: Members: jacint@140: jacint@140: T maxFlow() : returns the value of a maximum flow jacint@140: jacint@140: T flow(G_EdgeIt e) : for a fixed maximum flow x it returns x(e) jacint@140: jacint@140: G_EdgeMap_T_2 flow() : returns the fixed maximum flow x jacint@140: jacint@140: void minMinCut(G_NodeMap_bool& M) : jacint@140: sets M to the characteristic vector of the jacint@140: minimum min cut. M must be initialized to false. jacint@140: jacint@140: void maxMinCut(G_NodeMap_bool& M) : jacint@140: sets M to the characteristic vector of the jacint@140: maximum min cut. M must be initialized to false. jacint@140: jacint@140: void minCut(G_NodeMap_bool& M) : jacint@140: sets M to the characteristic vector of a jacint@140: min cut. M must be initialized to false. jacint@140: */ jacint@140: jacint@140: #ifndef PREFLOW_H jacint@140: #define PREFLOW_H jacint@140: jacint@140: #define H0 20 jacint@140: #define H1 1 jacint@140: jacint@140: #include jacint@140: #include jacint@140: jacint@140: #include //for error handling jacint@140: jacint@140: #include jacint@140: jacint@140: namespace hugo { jacint@140: jacint@140: template , jacint@140: typename FlowMap=typename Graph::EdgeMap > jacint@140: class PreflowAug { jacint@140: jacint@140: typedef typename Graph::NodeIt NodeIt; jacint@140: typedef typename Graph::EdgeIt EdgeIt; jacint@140: typedef typename Graph::EachNodeIt EachNodeIt; jacint@140: typedef typename Graph::EachEdgeIt EachEdgeIt; jacint@140: typedef typename Graph::OutEdgeIt OutEdgeIt; jacint@140: typedef typename Graph::InEdgeIt InEdgeIt; jacint@140: jacint@140: Graph& G; jacint@140: NodeIt s; jacint@140: NodeIt t; jacint@140: CapMap& capacity; jacint@140: FlowMap& _flow; jacint@140: bool flow_type; jacint@140: T value; jacint@140: jacint@140: public: jacint@140: double time; jacint@140: PreflowAug(Graph& _G, NodeIt _s, NodeIt _t, jacint@140: CapMap& _capacity, FlowMap& __flow, bool _flow_type ) : jacint@140: G(_G), s(_s), t(_t), capacity(_capacity), _flow(__flow), jacint@140: flow_type(_flow_type) jacint@140: { jacint@140: jacint@140: jacint@140: bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd jacint@140: int n=G.nodeNum(); jacint@140: int heur0=(int)(H0*n); //time while running 'bound decrease' jacint@140: int heur1=(int)(H1*n); //time while running 'highest label' jacint@140: int heur=heur1; //starting time interval (#of relabels) jacint@140: bool what_heur=1; jacint@140: /* jacint@140: what_heur is 0 in case 'bound decrease' jacint@140: and 1 in case 'highest label' jacint@140: */ jacint@140: bool end=false; jacint@140: /* jacint@140: Needed for 'bound decrease', 'true' jacint@140: means no active nodes are above bound b. jacint@140: */ jacint@140: int relabel=0; jacint@140: int k=n-2; //bound on the highest level under n containing a node jacint@140: int b=k; //bound on the highest level under n of an active node jacint@140: jacint@140: typename Graph::NodeMap level(G,n); jacint@140: typename Graph::NodeMap excess(G); jacint@140: jacint@140: std::vector active(n); jacint@140: typename Graph::NodeMap next(G); jacint@140: //Stack of the active nodes in level i < n. jacint@140: //We use it in both phases. jacint@140: jacint@140: typename Graph::NodeMap left(G); jacint@140: typename Graph::NodeMap right(G); jacint@140: std::vector level_list(n); jacint@140: /* jacint@140: List of the nodes in level i bfs_queue; jacint@140: bfs_queue.push(t); jacint@140: jacint@140: jacint@140: while (!bfs_queue.empty()) { jacint@140: jacint@140: NodeIt v=bfs_queue.front(); jacint@140: bfs_queue.pop(); jacint@140: int l=level.get(v)+1; jacint@140: jacint@140: for(InEdgeIt e=G.template first(v); e.valid(); ++e) { jacint@140: if ( capacity.get(e) == _flow.get(e) ) continue; jacint@140: NodeIt w=G.tail(e); jacint@140: if ( level.get(w) == n && w != s ) { jacint@140: bfs_queue.push(w); jacint@140: NodeIt first=level_list[l]; jacint@140: if ( first.valid() ) left.set(first,w); jacint@140: right.set(w,first); jacint@140: level_list[l]=w; jacint@140: level.set(w, l); jacint@140: } jacint@140: } jacint@140: jacint@140: for(OutEdgeIt e=G.template first(v); e.valid(); ++e) { jacint@140: if ( 0 == _flow.get(e) ) continue; jacint@140: NodeIt w=G.head(e); jacint@140: if ( level.get(w) == n && w != s ) { jacint@140: bfs_queue.push(w); jacint@140: NodeIt first=level_list[l]; jacint@140: if ( first != 0 ) left.set(first,w); jacint@140: right.set(w,first); jacint@140: level_list[l]=w; jacint@140: level.set(w, l); jacint@140: } jacint@140: } jacint@140: } jacint@140: jacint@140: level.set(s,n); jacint@140: jacint@140: jacint@140: jacint@140: /*The starting excess.*/ jacint@140: if ( flow_type ) { jacint@140: for(EachEdgeIt e=G.template first(); e.valid(); ++e) { jacint@140: if ( _flow.get(e) > 0 ) { jacint@140: T flo=_flow.get(e); jacint@140: NodeIt u=G.tail(e); jacint@140: NodeIt v=G.head(e); jacint@140: excess.set(u, excess.get(u)-flo); jacint@140: excess.set(v, excess.get(v)+flo); jacint@140: } jacint@140: } jacint@140: jacint@140: for(EachNodeIt v=G.template first(); v.valid(); ++v) { jacint@140: if ( excess.get(v) < 0 ) { jacint@140: std::cerr<<"It is not a pre_flow."<(s); e.valid(); ++e) jacint@140: { jacint@140: T c=capacity.get(e)-_flow.get(e); jacint@140: if ( c == 0 ) continue; jacint@140: NodeIt w=G.head(e); jacint@140: if ( level.get(w) < n ) { jacint@140: if ( excess.get(w) == 0 && w!=t ) { jacint@140: next.set(w,active[level.get(w)]); jacint@140: active[level.get(w)]=w; jacint@140: } jacint@140: _flow.set(e, capacity.get(e)); jacint@140: excess.set(w, excess.get(w)+c); jacint@140: } jacint@140: } jacint@140: jacint@140: /* jacint@140: End of preprocessing jacint@140: */ jacint@140: jacint@140: jacint@140: jacint@140: /* jacint@140: Push/relabel on the highest level active nodes. jacint@140: */ jacint@140: while ( true ) { jacint@140: jacint@140: if ( b == 0 ) { jacint@140: if ( phase ) break; jacint@140: jacint@140: if ( !what_heur && !end && k > 0 ) { jacint@140: b=k; jacint@140: end=true; jacint@140: } else { jacint@140: phase=1; jacint@140: time=currTime(); jacint@140: level.set(s,0); jacint@140: std::queue bfs_queue; jacint@140: bfs_queue.push(s); jacint@140: jacint@140: while (!bfs_queue.empty()) { jacint@140: jacint@140: NodeIt v=bfs_queue.front(); jacint@140: bfs_queue.pop(); jacint@140: int l=level.get(v)+1; jacint@140: jacint@140: for(InEdgeIt e=G.template first(v); e.valid(); ++e) { jacint@140: if ( capacity.get(e) == _flow.get(e) ) continue; jacint@140: NodeIt u=G.tail(e); jacint@140: if ( level.get(u) >= n ) { jacint@140: bfs_queue.push(u); jacint@140: level.set(u, l); jacint@140: if ( excess.get(u) > 0 ) { jacint@140: next.set(u,active[l]); jacint@140: active[l]=u; jacint@140: } jacint@140: } jacint@140: } jacint@140: jacint@140: for(OutEdgeIt e=G.template first(v); e.valid(); ++e) { jacint@140: if ( 0 == _flow.get(e) ) continue; jacint@140: NodeIt u=G.head(e); jacint@140: if ( level.get(u) >= n ) { jacint@140: bfs_queue.push(u); jacint@140: level.set(u, l); jacint@140: if ( excess.get(u) > 0 ) { jacint@140: next.set(u,active[l]); jacint@140: active[l]=u; jacint@140: } jacint@140: } jacint@140: } jacint@140: } jacint@140: b=n-2; jacint@140: } jacint@140: jacint@140: } jacint@140: jacint@140: jacint@140: if ( active[b] == 0 ) --b; jacint@140: else { jacint@140: end=false; jacint@140: jacint@140: NodeIt w=active[b]; jacint@140: active[b]=next.get(w); jacint@140: int lev=level.get(w); jacint@140: T exc=excess.get(w); jacint@140: int newlevel=n; //bound on the next level of w jacint@140: jacint@140: for(OutEdgeIt e=G.template first(w); e.valid(); ++e) { jacint@140: jacint@140: if ( _flow.get(e) == capacity.get(e) ) continue; jacint@140: NodeIt v=G.head(e); jacint@140: //e=wv jacint@140: jacint@140: if( lev > level.get(v) ) { jacint@140: /*Push is allowed now*/ jacint@140: jacint@140: if ( excess.get(v)==0 && v!=t && v!=s ) { jacint@140: int lev_v=level.get(v); jacint@140: next.set(v,active[lev_v]); jacint@140: active[lev_v]=v; jacint@140: } jacint@140: jacint@140: T cap=capacity.get(e); jacint@140: T flo=_flow.get(e); jacint@140: T remcap=cap-flo; jacint@140: jacint@140: if ( remcap >= exc ) { jacint@140: /*A nonsaturating push.*/ jacint@140: jacint@140: _flow.set(e, flo+exc); jacint@140: excess.set(v, excess.get(v)+exc); jacint@140: exc=0; jacint@140: break; jacint@140: jacint@140: } else { jacint@140: /*A saturating push.*/ jacint@140: jacint@140: _flow.set(e, cap); jacint@140: excess.set(v, excess.get(v)+remcap); jacint@140: exc-=remcap; jacint@140: } jacint@140: } else if ( newlevel > level.get(v) ){ jacint@140: newlevel = level.get(v); jacint@140: } jacint@140: jacint@140: } //for out edges wv jacint@140: jacint@140: jacint@140: if ( exc > 0 ) { jacint@140: for( InEdgeIt e=G.template first(w); e.valid(); ++e) { jacint@140: jacint@140: if( _flow.get(e) == 0 ) continue; jacint@140: NodeIt v=G.tail(e); jacint@140: //e=vw jacint@140: jacint@140: if( lev > level.get(v) ) { jacint@140: /*Push is allowed now*/ jacint@140: jacint@140: if ( excess.get(v)==0 && v!=t && v!=s ) { jacint@140: int lev_v=level.get(v); jacint@140: next.set(v,active[lev_v]); jacint@140: active[lev_v]=v; jacint@140: } jacint@140: jacint@140: T flo=_flow.get(e); jacint@140: jacint@140: if ( flo >= exc ) { jacint@140: /*A nonsaturating push.*/ jacint@140: jacint@140: _flow.set(e, flo-exc); jacint@140: excess.set(v, excess.get(v)+exc); jacint@140: exc=0; jacint@140: break; jacint@140: } else { jacint@140: /*A saturating push.*/ jacint@140: jacint@140: excess.set(v, excess.get(v)+flo); jacint@140: exc-=flo; jacint@140: _flow.set(e,0); jacint@140: } jacint@140: } else if ( newlevel > level.get(v) ) { jacint@140: newlevel = level.get(v); jacint@140: } jacint@140: } //for in edges vw jacint@140: jacint@140: } // if w still has excess after the out edge for cycle jacint@140: jacint@140: excess.set(w, exc); jacint@140: jacint@140: /* jacint@140: Relabel jacint@140: */ jacint@140: jacint@140: jacint@140: if ( exc > 0 ) { jacint@140: //now 'lev' is the old level of w jacint@140: jacint@140: if ( phase ) { jacint@140: level.set(w,++newlevel); jacint@140: next.set(w,active[newlevel]); jacint@140: active[newlevel]=w; jacint@140: b=newlevel; jacint@140: } else { jacint@140: //unlacing starts jacint@140: NodeIt right_n=right.get(w); jacint@140: NodeIt left_n=left.get(w); jacint@140: jacint@140: if ( right_n != 0 ) { jacint@140: if ( left_n != 0 ) { jacint@140: right.set(left_n, right_n); jacint@140: left.set(right_n, left_n); jacint@140: } else { jacint@140: level_list[lev]=right_n; jacint@140: left.set(right_n, 0); jacint@140: } jacint@140: } else { jacint@140: if ( left_n != 0 ) { jacint@140: right.set(left_n, 0); jacint@140: } else { jacint@140: level_list[lev]=0; jacint@140: jacint@140: } jacint@140: } jacint@140: //unlacing ends jacint@140: jacint@140: //gapping starts jacint@140: if ( level_list[lev]==0 ) { jacint@140: jacint@140: for (int i=lev; i!=k ; ) { jacint@140: NodeIt v=level_list[++i]; jacint@140: while ( v != 0 ) { jacint@140: level.set(v,n); jacint@140: v=right.get(v); jacint@140: } jacint@140: level_list[i]=0; jacint@140: if ( !what_heur ) active[i]=0; jacint@140: } jacint@140: jacint@140: level.set(w,n); jacint@140: b=lev-1; jacint@140: k=b; jacint@140: //gapping ends jacint@140: } else { jacint@140: jacint@140: if ( newlevel == n ) level.set(w,n); jacint@140: else { jacint@140: level.set(w,++newlevel); jacint@140: next.set(w,active[newlevel]); jacint@140: active[newlevel]=w; jacint@140: if ( what_heur ) b=newlevel; jacint@140: if ( k < newlevel ) ++k; jacint@140: NodeIt first=level_list[newlevel]; jacint@140: if ( first != 0 ) left.set(first,w); jacint@140: right.set(w,first); jacint@140: left.set(w,0); jacint@140: level_list[newlevel]=w; jacint@140: } jacint@140: } jacint@140: jacint@140: jacint@140: ++relabel; jacint@140: if ( relabel >= heur ) { jacint@140: relabel=0; jacint@140: if ( what_heur ) { jacint@140: what_heur=0; jacint@140: heur=heur0; jacint@140: end=false; jacint@140: } else { jacint@140: what_heur=1; jacint@140: heur=heur1; jacint@140: b=k; jacint@140: } jacint@140: } jacint@140: } //phase 0 jacint@140: jacint@140: jacint@140: } // if ( exc > 0 ) jacint@140: jacint@140: jacint@140: } // if stack[b] is nonempty jacint@140: jacint@140: } // while(true) jacint@140: jacint@140: jacint@140: value = excess.get(t); jacint@140: /*Max flow value.*/ jacint@140: jacint@140: } //void run() jacint@140: jacint@140: jacint@140: jacint@140: jacint@140: jacint@140: /* jacint@140: Returns the maximum value of a flow. jacint@140: */ jacint@140: jacint@140: T maxFlow() { jacint@140: return value; jacint@140: } jacint@140: jacint@140: jacint@140: jacint@140: /* jacint@140: For the maximum flow x found by the algorithm, jacint@140: it returns the flow value on edge e, i.e. x(e). jacint@140: */ jacint@140: jacint@140: T flow(EdgeIt e) { jacint@140: return _flow.get(e); jacint@140: } jacint@140: jacint@140: jacint@140: jacint@140: FlowMap flow() { jacint@140: return _flow; jacint@140: } jacint@140: jacint@140: jacint@140: jacint@140: void flow(FlowMap& __flow ) { jacint@140: for(EachNodeIt v=G.template first(); v.valid(); ++v) jacint@140: __flow.set(v,_flow.get(v)); jacint@140: } jacint@140: jacint@140: jacint@140: jacint@140: /* jacint@140: Returns the minimum min cut, by a bfs from s in the residual graph. jacint@140: */ jacint@140: jacint@140: template jacint@140: void minMinCut(_CutMap& M) { jacint@140: jacint@140: std::queue queue; jacint@140: jacint@140: M.set(s,true); jacint@140: queue.push(s); jacint@140: jacint@140: while (!queue.empty()) { jacint@140: NodeIt w=queue.front(); jacint@140: queue.pop(); jacint@140: jacint@140: for(OutEdgeIt e=G.template first(w) ; e.valid(); ++e) { jacint@140: NodeIt v=G.head(e); jacint@140: if (!M.get(v) && _flow.get(e) < capacity.get(e) ) { jacint@140: queue.push(v); jacint@140: M.set(v, true); jacint@140: } jacint@140: } jacint@140: jacint@140: for(InEdgeIt e=G.template first(w) ; e.valid(); ++e) { jacint@140: NodeIt v=G.tail(e); jacint@140: if (!M.get(v) && _flow.get(e) > 0 ) { jacint@140: queue.push(v); jacint@140: M.set(v, true); jacint@140: } jacint@140: } jacint@140: } jacint@140: } jacint@140: jacint@140: jacint@140: jacint@140: /* jacint@140: Returns the maximum min cut, by a reverse bfs jacint@140: from t in the residual graph. jacint@140: */ jacint@140: jacint@140: template jacint@140: void maxMinCut(_CutMap& M) { jacint@140: jacint@140: std::queue queue; jacint@140: jacint@140: M.set(t,true); jacint@140: queue.push(t); jacint@140: jacint@140: while (!queue.empty()) { jacint@140: NodeIt w=queue.front(); jacint@140: queue.pop(); jacint@140: jacint@140: for(InEdgeIt e=G.template first(w) ; e.valid(); ++e) { jacint@140: NodeIt v=G.tail(e); jacint@140: if (!M.get(v) && _flow.get(e) < capacity.get(e) ) { jacint@140: queue.push(v); jacint@140: M.set(v, true); jacint@140: } jacint@140: } jacint@140: jacint@140: for(OutEdgeIt e=G.template first(w) ; e.valid(); ++e) { jacint@140: NodeIt v=G.head(e); jacint@140: if (!M.get(v) && _flow.get(e) > 0 ) { jacint@140: queue.push(v); jacint@140: M.set(v, true); jacint@140: } jacint@140: } jacint@140: } jacint@140: jacint@140: for(EachNodeIt v=G.template first() ; v.valid(); ++v) { jacint@140: M.set(v, !M.get(v)); jacint@140: } jacint@140: jacint@140: } jacint@140: jacint@140: jacint@140: jacint@140: template jacint@140: void minCut(CutMap& M) { jacint@140: minMinCut(M); jacint@140: } jacint@140: jacint@140: jacint@140: }; jacint@140: }//namespace marci jacint@140: #endif jacint@140: jacint@140: jacint@140: jacint@140: