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