jacint@437: // -*- C++ -*- jacint@437: jacint@437: //run gyorsan tudna adni a minmincutot a 2 fazis elejen , ne vegyuk be konstruktorba egy cutmapet? jacint@437: //constzero jo igy? jacint@437: jacint@437: //majd marci megmondja betegyem-e bfs-t meg resgraphot jacint@437: jacint@437: //constzero helyett az kell hogy flow-e vagy csak preflow, ha flow akor csak jacint@437: //excess[t]-t kell szmaolni jacint@437: jacint@437: /* jacint@437: Heuristics: jacint@437: 2 phase jacint@437: gap jacint@437: list 'level_list' on the nodes on level i implemented by hand jacint@437: stack 'active' on the active nodes on level i implemented by hand jacint@437: runs heuristic 'highest label' for H1*n relabels jacint@437: runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label' jacint@437: jacint@437: Parameters H0 and H1 are initialized to 20 and 10. jacint@437: jacint@437: Constructors: jacint@437: jacint@437: Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if jacint@437: FlowMap is not constant zero, and should be true if it is jacint@437: jacint@437: Members: jacint@437: jacint@437: void run() jacint@437: jacint@437: T flowValue() : returns the value of a maximum flow jacint@437: jacint@437: void minMinCut(CutMap& M) : sets M to the characteristic vector of the jacint@437: minimum min cut. M should be a map of bools initialized to false. jacint@437: jacint@437: void maxMinCut(CutMap& M) : sets M to the characteristic vector of the jacint@437: maximum min cut. M should be a map of bools initialized to false. jacint@437: jacint@437: void minCut(CutMap& M) : sets M to the characteristic vector of jacint@437: a min cut. M should be a map of bools initialized to false. jacint@437: jacint@437: FIXME reset jacint@437: jacint@437: */ jacint@437: alpar@921: #ifndef LEMON_PREFLOW_H alpar@921: #define LEMON_PREFLOW_H jacint@437: jacint@437: #define H0 20 jacint@437: #define H1 1 jacint@437: jacint@437: #include jacint@437: #include jacint@437: #include jacint@437: alpar@921: namespace lemon { jacint@437: jacint@437: template , jacint@437: typename FlowMap=typename Graph::template EdgeMap > jacint@437: class Preflow { jacint@437: jacint@437: typedef typename Graph::Node Node; jacint@437: typedef typename Graph::Edge Edge; jacint@437: typedef typename Graph::NodeIt NodeIt; jacint@437: typedef typename Graph::OutEdgeIt OutEdgeIt; jacint@437: typedef typename Graph::InEdgeIt InEdgeIt; jacint@437: jacint@437: const Graph& G; jacint@437: Node s; jacint@437: Node t; jacint@437: const CapMap& capacity; jacint@437: FlowMap& flow; jacint@437: T value; jacint@437: bool constzero; jacint@437: bool isflow; jacint@437: jacint@437: public: jacint@437: Preflow(Graph& _G, Node _s, Node _t, CapMap& _capacity, jacint@437: FlowMap& _flow, bool _constzero, bool _isflow ) : jacint@437: G(_G), s(_s), t(_t), capacity(_capacity), flow(_flow), constzero(_constzero), isflow(_isflow) {} jacint@437: jacint@437: jacint@437: void run() { jacint@437: jacint@437: value=0; //for the subsequent runs jacint@437: jacint@437: bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd jacint@437: int n=G.nodeNum(); jacint@437: int heur0=(int)(H0*n); //time while running 'bound decrease' jacint@437: int heur1=(int)(H1*n); //time while running 'highest label' jacint@437: int heur=heur1; //starting time interval (#of relabels) jacint@437: bool what_heur=1; jacint@437: /* jacint@437: what_heur is 0 in case 'bound decrease' jacint@437: and 1 in case 'highest label' jacint@437: */ jacint@437: bool end=false; jacint@437: /* jacint@437: Needed for 'bound decrease', 'true' jacint@437: means no active nodes are above bound b. jacint@437: */ jacint@437: int relabel=0; jacint@437: int k=n-2; //bound on the highest level under n containing a node jacint@437: int b=k; //bound on the highest level under n of an active node jacint@437: jacint@437: typename Graph::template NodeMap level(G,n); jacint@437: typename Graph::template NodeMap excess(G); jacint@437: jacint@437: std::vector > active(n); jacint@437: /* std::vector active(n-1,INVALID); jacint@437: typename Graph::template NodeMap next(G,INVALID); jacint@437: //Stack of the active nodes in level i < n. jacint@437: //We use it in both phases.*/ jacint@437: jacint@437: typename Graph::template NodeMap left(G,INVALID); jacint@437: typename Graph::template NodeMap right(G,INVALID); jacint@437: std::vector level_list(n,INVALID); jacint@437: /* jacint@437: List of the nodes in level i bfs_queue; jacint@437: bfs_queue.push(t); jacint@437: jacint@437: while (!bfs_queue.empty()) { jacint@437: jacint@437: Node v=bfs_queue.front(); jacint@437: bfs_queue.pop(); jacint@437: int l=level[v]+1; jacint@437: jacint@437: InEdgeIt e; jacint@437: for(G.first(e,v); G.valid(e); G.next(e)) { jacint@437: Node w=G.tail(e); jacint@437: if ( level[w] == n && w != s ) { jacint@437: bfs_queue.push(w); jacint@437: Node first=level_list[l]; jacint@437: if ( G.valid(first) ) left.set(first,w); jacint@437: right.set(w,first); jacint@437: level_list[l]=w; jacint@437: level.set(w, l); jacint@437: } jacint@437: } jacint@437: } jacint@437: jacint@437: //the starting flow jacint@437: OutEdgeIt e; jacint@437: for(G.first(e,s); G.valid(e); G.next(e)) jacint@437: { jacint@437: T c=capacity[e]; jacint@437: if ( c == 0 ) continue; jacint@437: Node w=G.head(e); jacint@437: if ( level[w] < n ) { jacint@437: if ( excess[w] == 0 && w!=t ) active[level[w]].push(w); jacint@437: flow.set(e, c); jacint@437: excess.set(w, excess[w]+c); jacint@437: } jacint@437: } jacint@437: } jacint@437: else jacint@437: { jacint@437: jacint@437: /* jacint@437: Reverse_bfs from t in the residual graph, jacint@437: to find the starting level. jacint@437: */ jacint@437: level.set(t,0); jacint@437: std::queue bfs_queue; jacint@437: bfs_queue.push(t); jacint@437: jacint@437: while (!bfs_queue.empty()) { jacint@437: jacint@437: Node v=bfs_queue.front(); jacint@437: bfs_queue.pop(); jacint@437: int l=level[v]+1; jacint@437: jacint@437: InEdgeIt e; jacint@437: for(G.first(e,v); G.valid(e); G.next(e)) { jacint@437: if ( capacity[e] == flow[e] ) continue; jacint@437: Node w=G.tail(e); jacint@437: if ( level[w] == n && w != s ) { jacint@437: bfs_queue.push(w); jacint@437: Node first=level_list[l]; jacint@437: if ( G.valid(first) ) left.set(first,w); jacint@437: right.set(w,first); jacint@437: level_list[l]=w; jacint@437: level.set(w, l); jacint@437: } jacint@437: } jacint@437: jacint@437: OutEdgeIt f; jacint@437: for(G.first(f,v); G.valid(f); G.next(f)) { jacint@437: if ( 0 == flow[f] ) continue; jacint@437: Node w=G.head(f); jacint@437: if ( level[w] == n && w != s ) { jacint@437: bfs_queue.push(w); jacint@437: Node first=level_list[l]; jacint@437: if ( G.valid(first) ) left.set(first,w); jacint@437: right.set(w,first); jacint@437: level_list[l]=w; jacint@437: level.set(w, l); jacint@437: } jacint@437: } jacint@437: } jacint@437: jacint@437: jacint@437: /* jacint@437: Counting the excess jacint@437: */ jacint@437: jacint@437: if ( !isflow ) { jacint@437: NodeIt v; jacint@437: for(G.first(v); G.valid(v); G.next(v)) { jacint@437: T exc=0; jacint@437: jacint@437: InEdgeIt e; jacint@437: for(G.first(e,v); G.valid(e); G.next(e)) exc+=flow[e]; jacint@437: OutEdgeIt f; jacint@437: for(G.first(f,v); G.valid(f); G.next(f)) exc-=flow[f]; jacint@437: jacint@437: excess.set(v,exc); jacint@437: jacint@437: //putting the active nodes into the stack jacint@437: int lev=level[v]; jacint@437: if ( exc > 0 && lev < n && v != t ) active[lev].push(v); jacint@437: } jacint@437: } else { jacint@437: T exc=0; jacint@437: jacint@437: InEdgeIt e; jacint@437: for(G.first(e,t); G.valid(e); G.next(e)) exc+=flow[e]; jacint@437: OutEdgeIt f; jacint@437: for(G.first(f,t); G.valid(f); G.next(f)) exc-=flow[f]; jacint@437: jacint@437: excess.set(t,exc); jacint@437: } jacint@437: jacint@437: jacint@437: //the starting flow jacint@437: OutEdgeIt e; jacint@437: for(G.first(e,s); G.valid(e); G.next(e)) jacint@437: { jacint@437: T rem=capacity[e]-flow[e]; jacint@437: if ( rem == 0 ) continue; jacint@437: Node w=G.head(e); jacint@437: if ( level[w] < n ) { jacint@437: if ( excess[w] == 0 && w!=t ) active[level[w]].push(w); jacint@437: flow.set(e, capacity[e]); jacint@437: excess.set(w, excess[w]+rem); jacint@437: } jacint@437: } jacint@437: jacint@437: InEdgeIt f; jacint@437: for(G.first(f,s); G.valid(f); G.next(f)) jacint@437: { jacint@437: if ( flow[f] == 0 ) continue; jacint@437: Node w=G.tail(f); jacint@437: if ( level[w] < n ) { jacint@437: if ( excess[w] == 0 && w!=t ) active[level[w]].push(w); jacint@437: excess.set(w, excess[w]+flow[f]); jacint@437: flow.set(f, 0); jacint@437: } jacint@437: } jacint@437: } jacint@437: jacint@437: jacint@437: jacint@437: jacint@437: /* jacint@437: End of preprocessing jacint@437: */ jacint@437: jacint@437: jacint@437: jacint@437: /* jacint@437: Push/relabel on the highest level active nodes. jacint@437: */ jacint@437: while ( true ) { jacint@437: jacint@437: if ( b == 0 ) { jacint@437: if ( phase ) break; jacint@437: jacint@437: if ( !what_heur && !end && k > 0 ) { jacint@437: b=k; jacint@437: end=true; jacint@437: } else { jacint@437: phase=1; jacint@437: level.set(s,0); jacint@437: std::queue bfs_queue; jacint@437: bfs_queue.push(s); jacint@437: jacint@437: while (!bfs_queue.empty()) { jacint@437: jacint@437: Node v=bfs_queue.front(); jacint@437: bfs_queue.pop(); jacint@437: int l=level[v]+1; jacint@437: jacint@437: InEdgeIt e; jacint@437: for(G.first(e,v); G.valid(e); G.next(e)) { jacint@437: if ( capacity[e] == flow[e] ) continue; jacint@437: Node u=G.tail(e); jacint@437: if ( level[u] >= n ) { jacint@437: bfs_queue.push(u); jacint@437: level.set(u, l); jacint@437: if ( excess[u] > 0 ) active[l].push(u); jacint@437: } jacint@437: } jacint@437: jacint@437: OutEdgeIt f; jacint@437: for(G.first(f,v); G.valid(f); G.next(f)) { jacint@437: if ( 0 == flow[f] ) continue; jacint@437: Node u=G.head(f); jacint@437: if ( level[u] >= n ) { jacint@437: bfs_queue.push(u); jacint@437: level.set(u, l); jacint@437: if ( excess[u] > 0 ) active[l].push(u); jacint@437: } jacint@437: } jacint@437: } jacint@437: b=n-2; jacint@437: } jacint@437: jacint@437: } jacint@437: jacint@437: jacint@437: /// jacint@437: if ( active[b].empty() ) --b; jacint@437: else { jacint@437: end=false; jacint@437: jacint@437: Node w=active[b].top(); jacint@437: active[b].pop(); jacint@437: int lev=level[w]; jacint@437: T exc=excess[w]; jacint@437: int newlevel=n; //bound on the next level of w jacint@437: jacint@437: OutEdgeIt e; jacint@437: for(G.first(e,w); G.valid(e); G.next(e)) { jacint@437: jacint@437: if ( flow[e] == capacity[e] ) continue; jacint@437: Node v=G.head(e); jacint@437: //e=wv jacint@437: jacint@437: if( lev > level[v] ) { jacint@437: /*Push is allowed now*/ jacint@437: jacint@437: if ( excess[v]==0 && v!=t && v!=s ) { jacint@437: int lev_v=level[v]; jacint@437: active[lev_v].push(v); jacint@437: } jacint@437: jacint@437: T cap=capacity[e]; jacint@437: T flo=flow[e]; jacint@437: T remcap=cap-flo; jacint@437: jacint@437: if ( remcap >= exc ) { jacint@437: /*A nonsaturating push.*/ jacint@437: jacint@437: flow.set(e, flo+exc); jacint@437: excess.set(v, excess[v]+exc); jacint@437: exc=0; jacint@437: break; jacint@437: jacint@437: } else { jacint@437: /*A saturating push.*/ jacint@437: jacint@437: flow.set(e, cap); jacint@437: excess.set(v, excess[v]+remcap); jacint@437: exc-=remcap; jacint@437: } jacint@437: } else if ( newlevel > level[v] ){ jacint@437: newlevel = level[v]; jacint@437: } jacint@437: jacint@437: } //for out edges wv jacint@437: jacint@437: jacint@437: if ( exc > 0 ) { jacint@437: InEdgeIt e; jacint@437: for(G.first(e,w); G.valid(e); G.next(e)) { jacint@437: jacint@437: if( flow[e] == 0 ) continue; jacint@437: Node v=G.tail(e); jacint@437: //e=vw jacint@437: jacint@437: if( lev > level[v] ) { jacint@437: /*Push is allowed now*/ jacint@437: jacint@437: if ( excess[v]==0 && v!=t && v!=s ) { jacint@437: int lev_v=level[v]; jacint@437: active[lev_v].push(v); jacint@437: } jacint@437: jacint@437: T flo=flow[e]; jacint@437: jacint@437: if ( flo >= exc ) { jacint@437: /*A nonsaturating push.*/ jacint@437: jacint@437: flow.set(e, flo-exc); jacint@437: excess.set(v, excess[v]+exc); jacint@437: exc=0; jacint@437: break; jacint@437: } else { jacint@437: /*A saturating push.*/ jacint@437: jacint@437: excess.set(v, excess[v]+flo); jacint@437: exc-=flo; jacint@437: flow.set(e,0); jacint@437: } jacint@437: } else if ( newlevel > level[v] ) { jacint@437: newlevel = level[v]; jacint@437: } jacint@437: } //for in edges vw jacint@437: jacint@437: } // if w still has excess after the out edge for cycle jacint@437: jacint@437: excess.set(w, exc); jacint@437: /// push jacint@437: jacint@437: jacint@437: /* jacint@437: Relabel jacint@437: */ jacint@437: jacint@437: jacint@437: if ( exc > 0 ) { jacint@437: //now 'lev' is the old level of w jacint@437: jacint@437: if ( phase ) { jacint@437: level.set(w,++newlevel); jacint@437: active[newlevel].push(w); jacint@437: b=newlevel; jacint@437: } else { jacint@437: //unlacing starts jacint@437: Node right_n=right[w]; jacint@437: Node left_n=left[w]; jacint@437: jacint@437: if ( G.valid(right_n) ) { jacint@437: if ( G.valid(left_n) ) { jacint@437: right.set(left_n, right_n); jacint@437: left.set(right_n, left_n); jacint@437: } else { jacint@437: level_list[lev]=right_n; jacint@437: left.set(right_n, INVALID); jacint@437: } jacint@437: } else { jacint@437: if ( G.valid(left_n) ) { jacint@437: right.set(left_n, INVALID); jacint@437: } else { jacint@437: level_list[lev]=INVALID; jacint@437: } jacint@437: } jacint@437: //unlacing ends jacint@437: jacint@437: if ( !G.valid(level_list[lev]) ) { jacint@437: jacint@437: //gapping starts jacint@437: for (int i=lev; i!=k ; ) { jacint@437: Node v=level_list[++i]; jacint@437: while ( G.valid(v) ) { jacint@437: level.set(v,n); jacint@437: v=right[v]; jacint@437: } jacint@437: level_list[i]=INVALID; jacint@437: if ( !what_heur ) { jacint@437: while ( !active[i].empty() ) { jacint@437: active[i].pop(); //FIXME: ezt szebben kene jacint@437: } jacint@437: } jacint@437: } jacint@437: jacint@437: level.set(w,n); jacint@437: b=lev-1; jacint@437: k=b; jacint@437: //gapping ends jacint@437: jacint@437: } else { jacint@437: jacint@437: if ( newlevel == n ) level.set(w,n); jacint@437: else { jacint@437: level.set(w,++newlevel); jacint@437: active[newlevel].push(w); jacint@437: if ( what_heur ) b=newlevel; jacint@437: if ( k < newlevel ) ++k; //now k=newlevel jacint@437: Node first=level_list[newlevel]; jacint@437: if ( G.valid(first) ) left.set(first,w); jacint@437: right.set(w,first); jacint@437: left.set(w,INVALID); jacint@437: level_list[newlevel]=w; jacint@437: } jacint@437: } jacint@437: jacint@437: jacint@437: ++relabel; jacint@437: if ( relabel >= heur ) { jacint@437: relabel=0; jacint@437: if ( what_heur ) { jacint@437: what_heur=0; jacint@437: heur=heur0; jacint@437: end=false; jacint@437: } else { jacint@437: what_heur=1; jacint@437: heur=heur1; jacint@437: b=k; jacint@437: } jacint@437: } jacint@437: } //phase 0 jacint@437: jacint@437: jacint@437: } // if ( exc > 0 ) jacint@437: jacint@437: jacint@437: } // if stack[b] is nonempty jacint@437: jacint@437: } // while(true) jacint@437: jacint@437: jacint@437: value = excess[t]; jacint@437: /*Max flow value.*/ jacint@437: jacint@437: } //void run() jacint@437: jacint@437: jacint@437: jacint@437: jacint@437: jacint@437: /* jacint@437: Returns the maximum value of a flow. jacint@437: */ jacint@437: jacint@437: T flowValue() { jacint@437: return value; jacint@437: } jacint@437: jacint@437: jacint@437: FlowMap Flow() { jacint@437: return flow; jacint@437: } jacint@437: jacint@437: jacint@437: void Flow(FlowMap& _flow ) { jacint@437: NodeIt v; jacint@437: for(G.first(v) ; G.valid(v); G.next(v)) jacint@437: _flow.set(v,flow[v]); jacint@437: } jacint@437: jacint@437: jacint@437: jacint@437: /* jacint@437: Returns the minimum min cut, by a bfs from s in the residual graph. jacint@437: */ jacint@437: jacint@437: template jacint@437: void minMinCut(_CutMap& M) { jacint@437: jacint@437: std::queue queue; jacint@437: jacint@437: M.set(s,true); jacint@437: queue.push(s); jacint@437: jacint@437: while (!queue.empty()) { jacint@437: Node w=queue.front(); jacint@437: queue.pop(); jacint@437: jacint@437: OutEdgeIt e; jacint@437: for(G.first(e,w) ; G.valid(e); G.next(e)) { jacint@437: Node v=G.head(e); jacint@437: if (!M[v] && flow[e] < capacity[e] ) { jacint@437: queue.push(v); jacint@437: M.set(v, true); jacint@437: } jacint@437: } jacint@437: jacint@437: InEdgeIt f; jacint@437: for(G.first(f,w) ; G.valid(f); G.next(f)) { jacint@437: Node v=G.tail(f); jacint@437: if (!M[v] && flow[f] > 0 ) { jacint@437: queue.push(v); jacint@437: M.set(v, true); jacint@437: } jacint@437: } jacint@437: } jacint@437: } jacint@437: jacint@437: jacint@437: jacint@437: /* jacint@437: Returns the maximum min cut, by a reverse bfs jacint@437: from t in the residual graph. jacint@437: */ jacint@437: jacint@437: template jacint@437: void maxMinCut(_CutMap& M) { jacint@437: jacint@437: std::queue queue; jacint@437: jacint@437: M.set(t,true); jacint@437: queue.push(t); jacint@437: jacint@437: while (!queue.empty()) { jacint@437: Node w=queue.front(); jacint@437: queue.pop(); jacint@437: jacint@437: jacint@437: InEdgeIt e; jacint@437: for(G.first(e,w) ; G.valid(e); G.next(e)) { jacint@437: Node v=G.tail(e); jacint@437: if (!M[v] && flow[e] < capacity[e] ) { jacint@437: queue.push(v); jacint@437: M.set(v, true); jacint@437: } jacint@437: } jacint@437: jacint@437: OutEdgeIt f; jacint@437: for(G.first(f,w) ; G.valid(f); G.next(f)) { jacint@437: Node v=G.head(f); jacint@437: if (!M[v] && flow[f] > 0 ) { jacint@437: queue.push(v); jacint@437: M.set(v, true); jacint@437: } jacint@437: } jacint@437: } jacint@437: jacint@437: NodeIt v; jacint@437: for(G.first(v) ; G.valid(v); G.next(v)) { jacint@437: M.set(v, !M[v]); jacint@437: } jacint@437: jacint@437: } jacint@437: jacint@437: jacint@437: jacint@437: template jacint@437: void minCut(CutMap& M) { jacint@437: minMinCut(M); jacint@437: } jacint@437: jacint@437: jacint@437: void resetTarget (Node _t) {t=_t;} jacint@437: void resetSource (Node _s) {s=_s;} jacint@437: jacint@437: void resetCap (CapMap _cap) {capacity=_cap;} jacint@437: jacint@437: void resetFlow (FlowMap _flow, bool _constzero) { jacint@437: flow=_flow; jacint@437: constzero=_constzero; jacint@437: } jacint@437: jacint@437: jacint@437: jacint@437: }; jacint@437: alpar@921: } //namespace lemon jacint@437: jacint@437: #endif //PREFLOW_H jacint@437: jacint@437: jacint@437: jacint@437: