// -*- C++ -*- /* preflow_push_hl.h by jacint. Runs the highest label variant of the preflow push algorithm with running time O(n^2\sqrt(m)). Member functions: void run() : runs the algorithm The following functions should be used after run() was already run. T maxflow() : returns the value of a maximum flow T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e) Graph::EdgeMap allflow() : returns the fixed maximum flow x Graph::NodeMap mincut() : returns a characteristic vector of a minimum cut. (An empty level in the algorithm gives a minimum cut.) */ #ifndef PREFLOW_PUSH_HL_H #define PREFLOW_PUSH_HL_H #define A 1 #include #include #include namespace marci { template class preflow_push_hl { typedef typename Graph::NodeIt NodeIt; typedef typename Graph::EdgeIt EdgeIt; typedef typename Graph::EachNodeIt EachNodeIt; typedef typename Graph::OutEdgeIt OutEdgeIt; typedef typename Graph::InEdgeIt InEdgeIt; Graph& G; NodeIt s; NodeIt t; typename Graph::EdgeMap flow; typename Graph::EdgeMap capacity; T value; typename Graph::NodeMap mincutvector; public: preflow_push_hl(Graph& _G, NodeIt _s, NodeIt _t, typename Graph::EdgeMap& _capacity) : G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity), mincutvector(_G, true) { } /* The run() function runs the highest label preflow-push, running time: O(n^2\sqrt(m)) */ void run() { std::cout<<"A is "< level(G); typename Graph::NodeMap excess(G); int n=G.nodeNum(); int b=n-2; /* b is a bound on the highest level of an active node. In the beginning it is at most n-2. */ std::vector numb(n); //The number of nodes on level i < n. std::vector > stack(2*n-1); //Stack of the active nodes in level i. /*Reverse_bfs from t, to find the starting level.*/ reverse_bfs bfs(G, t); bfs.run(); for(EachNodeIt v=G.template first(); v.valid(); ++v) { int dist=bfs.dist(v); level.set(v, dist); ++numb[dist]; } level.set(s,n); /* Starting flow. It is everywhere 0 at the moment. */ for(OutEdgeIt e=G.template first(s); e.valid(); ++e) { if ( capacity.get(e) > 0 ) { NodeIt w=G.head(e); if ( w!=s ) { if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w); flow.set(e, capacity.get(e)); excess.set(w, excess.get(w)+capacity.get(e)); } } } /* End of preprocessing */ /* Push/relabel on the highest level active nodes. */ /*While there exists an active node.*/ while (b) { /*We decrease the bound if there is no active node of level b.*/ if (stack[b].empty()) { --b; } else { NodeIt w=stack[b].top(); //w is a highest label active node. stack[b].pop(); int newlevel=2*n-2; //In newlevel we bound the next level of w. for(OutEdgeIt e=G.template first(w); e.valid(); ++e) { if ( flow.get(e) < capacity.get(e) ) { /*e is an edge of the residual graph */ NodeIt v=G.head(e); /*e is the edge wv.*/ if( level.get(w) == level.get(v)+1 ) { /*Push is allowed now*/ if ( excess.get(v)==0 && v != s && v !=t ) stack[level.get(v)].push(v); /*v becomes active.*/ if ( capacity.get(e)-flow.get(e) > excess.get(w) ) { /*A nonsaturating push.*/ flow.set(e, flow.get(e)+excess.get(w)); excess.set(v, excess.get(v)+excess.get(w)); excess.set(w,0); break; } else { /*A saturating push.*/ excess.set(v, excess.get(v)+capacity.get(e)-flow.get(e)); excess.set(w, excess.get(w)-capacity.get(e)+flow.get(e)); flow.set(e, capacity.get(e)); if ( excess.get(w)==0 ) break; /*If w is not active any more, then we go on to the next node.*/ } } else { newlevel = newlevel < level.get(v) ? newlevel : level.get(v); } } //if the out edge wv is in the res graph } //for out edges wv if ( excess.get(w) > 0 ) { for( InEdgeIt e=G.template first(w); e.valid(); ++e) { NodeIt v=G.tail(e); /*e is the edge vw.*/ if( flow.get(e) > 0 ) { /*e is an edge of the residual graph */ if( level.get(w)==level.get(v)+1 ) { /*Push is allowed now*/ if ( excess.get(v)==0 && v != s && v !=t) stack[level.get(v)].push(v); /*v becomes active.*/ if ( flow.get(e) > excess.get(w) ) { /*A nonsaturating push.*/ flow.set(e, flow.get(e)-excess.get(w)); excess.set(v, excess.get(v)+excess.get(w)); excess.set(w,0); break; } else { /*A saturating push.*/ excess.set(v, excess.get(v)+flow.get(e)); excess.set(w, excess.get(w)-flow.get(e)); flow.set(e,0); if ( excess.get(w)==0 ) break; } } else { newlevel = newlevel < level.get(v) ? newlevel : level.get(v); } } //if in edge vw is in the res graph } //for in edges vw } // if w still has excess after the out edge for cycle /* Relabel */ if ( excess.get(w) > 0 ) { int oldlevel=level.get(w); level.set(w,++newlevel); if ( oldlevel < n ) { --numb[oldlevel]; if ( !numb[oldlevel] && oldlevel < A*n ) { //If the level of w gets empty. for (EachNodeIt v=G.template first(); v.valid() ; ++v) { if (level.get(v) > oldlevel && level.get(v) < n ) level.set(v,n); } for (int i=oldlevel+1 ; i!=n ; ++i) numb[i]=0; if ( newlevel < n ) newlevel=n; } else { if ( newlevel < n ) ++numb[newlevel]; } } else { if ( newlevel < n ) ++numb[newlevel]; } stack[newlevel].push(w); b=newlevel; } } // if stack[b] is nonempty } // while(b) value = excess.get(t); /*Max flow value.*/ } //void run() /* Returns the maximum value of a flow. */ T maxflow() { return value; } /* For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e). */ T flowonedge(EdgeIt e) { return flow.get(e); } /* Returns the maximum flow x found by the algorithm. */ typename Graph::EdgeMap allflow() { return flow; } /* Returns a minimum cut by using a reverse bfs from t in the residual graph. */ typename Graph::NodeMap mincut() { std::queue queue; mincutvector.set(t,false); queue.push(t); while (!queue.empty()) { NodeIt w=queue.front(); queue.pop(); for(InEdgeIt e=G.template first(w) ; e.valid(); ++e) { NodeIt v=G.tail(e); if (mincutvector.get(v) && flow.get(e) < capacity.get(e) ) { queue.push(v); mincutvector.set(v, false); } } // for for(OutEdgeIt e=G.template first(w) ; e.valid(); ++e) { NodeIt v=G.head(e); if (mincutvector.get(v) && flow.get(e) > 0 ) { queue.push(v); mincutvector.set(v, false); } } // for } return mincutvector; } }; }//namespace marci #endif