1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/jacint/preflow_hl2.h Wed Feb 18 14:43:01 2004 +0000
1.3 @@ -0,0 +1,403 @@
1.4 +// -*- C++ -*-
1.5 +/*
1.6 +preflow_hl2.h
1.7 +by jacint.
1.8 +Runs the highest label variant of the preflow push algorithm with
1.9 +running time O(n^2\sqrt(m)), with the 'empty level' and with the
1.10 +heuristic that the bound b on the active nodes is not increased
1.11 +only when b=0, when we put b=2*n-2.
1.12 +
1.13 +'A' is a parameter for the empty_level heuristic
1.14 +
1.15 +Member functions:
1.16 +
1.17 +void run() : runs the algorithm
1.18 +
1.19 + The following functions should be used after run() was already run.
1.20 +
1.21 +T maxflow() : returns the value of a maximum flow
1.22 +
1.23 +T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
1.24 +
1.25 +FlowMap allflow() : returns the fixed maximum flow x
1.26 +
1.27 +void mincut(CutMap& M) : sets M to the characteristic vector of a
1.28 + minimum cut. M should be a map of bools initialized to false.
1.29 +
1.30 +void min_mincut(CutMap& M) : sets M to the characteristic vector of the
1.31 + minimum min cut. M should be a map of bools initialized to false.
1.32 +
1.33 +void max_mincut(CutMap& M) : sets M to the characteristic vector of the
1.34 + maximum min cut. M should be a map of bools initialized to false.
1.35 +
1.36 +*/
1.37 +
1.38 +#ifndef PREFLOW_HL2_H
1.39 +#define PREFLOW_HL2_H
1.40 +
1.41 +#define A 1
1.42 +
1.43 +#include <vector>
1.44 +#include <stack>
1.45 +#include <queue>
1.46 +
1.47 +namespace marci {
1.48 +
1.49 + template <typename Graph, typename T,
1.50 + typename FlowMap=typename Graph::EdgeMap<T>, typename CapMap=typename Graph::EdgeMap<T>,
1.51 + typename IntMap=typename Graph::NodeMap<int>, typename TMap=typename Graph::NodeMap<T> >
1.52 + class preflow_hl2 {
1.53 +
1.54 + typedef typename Graph::NodeIt NodeIt;
1.55 + typedef typename Graph::EdgeIt EdgeIt;
1.56 + typedef typename Graph::EachNodeIt EachNodeIt;
1.57 + typedef typename Graph::OutEdgeIt OutEdgeIt;
1.58 + typedef typename Graph::InEdgeIt InEdgeIt;
1.59 +
1.60 + Graph& G;
1.61 + NodeIt s;
1.62 + NodeIt t;
1.63 + FlowMap flow;
1.64 + CapMap& capacity;
1.65 + T value;
1.66 +
1.67 + public:
1.68 +
1.69 + preflow_hl2(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
1.70 + G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) { }
1.71 +
1.72 +
1.73 + void run() {
1.74 +
1.75 + bool no_end=true;
1.76 + int n=G.nodeNum();
1.77 + int b=n-2;
1.78 + /*
1.79 + b is a bound on the highest level of an active node.
1.80 + In the beginning it is at most n-2.
1.81 + */
1.82 +
1.83 + IntMap level(G,n);
1.84 + TMap excess(G);
1.85 +
1.86 + std::vector<int> numb(n+1);
1.87 + /*
1.88 + The number of nodes on level i < n. It is
1.89 + initialized to n+1, because of the reverse_bfs-part.
1.90 + */
1.91 +
1.92 + std::vector<std::stack<NodeIt> > stack(2*n-1);
1.93 + //Stack of the active nodes in level i.
1.94 +
1.95 +
1.96 + /*Reverse_bfs from t, to find the starting level.*/
1.97 + level.set(t,0);
1.98 + std::queue<NodeIt> bfs_queue;
1.99 + bfs_queue.push(t);
1.100 +
1.101 + while (!bfs_queue.empty()) {
1.102 +
1.103 + NodeIt v=bfs_queue.front();
1.104 + bfs_queue.pop();
1.105 + int l=level.get(v)+1;
1.106 +
1.107 + for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
1.108 + NodeIt w=G.tail(e);
1.109 + if ( level.get(w) == n ) {
1.110 + bfs_queue.push(w);
1.111 + ++numb[l];
1.112 + level.set(w, l);
1.113 + }
1.114 + }
1.115 + }
1.116 +
1.117 + level.set(s,n);
1.118 +
1.119 +
1.120 +
1.121 + /* Starting flow. It is everywhere 0 at the moment. */
1.122 + for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
1.123 + {
1.124 + if ( capacity.get(e) == 0 ) continue;
1.125 + NodeIt w=G.head(e);
1.126 + if ( w!=s ) {
1.127 + if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w);
1.128 + flow.set(e, capacity.get(e));
1.129 + excess.set(w, excess.get(w)+capacity.get(e));
1.130 + }
1.131 + }
1.132 +
1.133 + /*
1.134 + End of preprocessing
1.135 + */
1.136 +
1.137 +
1.138 +
1.139 + /*
1.140 + Push/relabel on the highest level active nodes.
1.141 + */
1.142 + /*While there exists an active node.*/
1.143 + while (b) {
1.144 + if ( stack[b].empty() ) {
1.145 + if ( b==1 ) {
1.146 + if ( !no_end ) break;
1.147 + else {
1.148 + b=2*n-2;
1.149 + no_end=false;
1.150 + }
1.151 + }
1.152 + --b;
1.153 + } else {
1.154 +
1.155 + no_end=true;
1.156 +
1.157 + NodeIt w=stack[b].top(); //w is a highest label active node.
1.158 + stack[b].pop();
1.159 + int lev=level.get(w);
1.160 + int exc=excess.get(w);
1.161 + int newlevel=2*n-2; //In newlevel we bound the next level of w.
1.162 +
1.163 + // if ( level.get(w) < n ) { //Nem tudom ez mukodik-e
1.164 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
1.165 +
1.166 + if ( flow.get(e) == capacity.get(e) ) continue;
1.167 + NodeIt v=G.head(e);
1.168 + //e=wv
1.169 +
1.170 + if( lev > level.get(v) ) {
1.171 + /*Push is allowed now*/
1.172 +
1.173 + if ( excess.get(v)==0 && v != s && v !=t )
1.174 + stack[level.get(v)].push(v);
1.175 + /*v becomes active.*/
1.176 +
1.177 + int cap=capacity.get(e);
1.178 + int flo=flow.get(e);
1.179 + int remcap=cap-flo;
1.180 +
1.181 + if ( remcap >= exc ) {
1.182 + /*A nonsaturating push.*/
1.183 +
1.184 + flow.set(e, flo+exc);
1.185 + excess.set(v, excess.get(v)+exc);
1.186 + exc=0;
1.187 + break;
1.188 +
1.189 + } else {
1.190 + /*A saturating push.*/
1.191 +
1.192 + flow.set(e, cap );
1.193 + excess.set(v, excess.get(v)+remcap);
1.194 + exc-=remcap;
1.195 + }
1.196 + } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
1.197 +
1.198 + } //for out edges wv
1.199 +
1.200 +
1.201 + if ( exc > 0 ) {
1.202 + for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
1.203 +
1.204 + if( flow.get(e) == 0 ) continue;
1.205 + NodeIt v=G.tail(e);
1.206 + //e=vw
1.207 +
1.208 + if( lev > level.get(v) ) {
1.209 + /*Push is allowed now*/
1.210 +
1.211 + if ( excess.get(v)==0 && v != s && v !=t)
1.212 + stack[level.get(v)].push(v);
1.213 + /*v becomes active.*/
1.214 +
1.215 + int flo=flow.get(e);
1.216 +
1.217 + if ( flo >= exc ) {
1.218 + /*A nonsaturating push.*/
1.219 +
1.220 + flow.set(e, flo-exc);
1.221 + excess.set(v, excess.get(v)+exc);
1.222 + exc=0;
1.223 + break;
1.224 + } else {
1.225 + /*A saturating push.*/
1.226 +
1.227 + excess.set(v, excess.get(v)+flo);
1.228 + exc-=flo;
1.229 + flow.set(e,0);
1.230 + }
1.231 + } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
1.232 +
1.233 + } //for in edges vw
1.234 +
1.235 + } // if w still has excess after the out edge for cycle
1.236 +
1.237 + excess.set(w, exc);
1.238 +
1.239 +
1.240 + /*
1.241 + Relabel
1.242 + */
1.243 +
1.244 + if ( exc > 0 ) {
1.245 + //now 'lev' is the old level of w
1.246 + level.set(w,++newlevel);
1.247 +
1.248 + if ( lev < n ) {
1.249 + --numb[lev];
1.250 +
1.251 + if ( !numb[lev] && lev < A*n ) { //If the level of w gets empty.
1.252 +
1.253 + for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
1.254 + if (level.get(v) > lev && level.get(v) < n ) level.set(v,n);
1.255 + }
1.256 + for (int i=lev+1 ; i!=n ; ++i) numb[i]=0;
1.257 + if ( newlevel < n ) newlevel=n;
1.258 + } else {
1.259 + if ( newlevel < n ) ++numb[newlevel];
1.260 + }
1.261 + }
1.262 +
1.263 + stack[newlevel].push(w);
1.264 +
1.265 + }
1.266 +
1.267 + } // if stack[b] is nonempty
1.268 +
1.269 + } // while(b)
1.270 +
1.271 +
1.272 + value = excess.get(t);
1.273 + /*Max flow value.*/
1.274 +
1.275 +
1.276 + } //void run()
1.277 +
1.278 +
1.279 +
1.280 +
1.281 +
1.282 + /*
1.283 + Returns the maximum value of a flow.
1.284 + */
1.285 +
1.286 + T maxflow() {
1.287 + return value;
1.288 + }
1.289 +
1.290 +
1.291 +
1.292 + /*
1.293 + For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e).
1.294 + */
1.295 +
1.296 + T flowonedge(EdgeIt e) {
1.297 + return flow.get(e);
1.298 + }
1.299 +
1.300 +
1.301 +
1.302 + /*
1.303 + Returns the maximum flow x found by the algorithm.
1.304 + */
1.305 +
1.306 + FlowMap allflow() {
1.307 + return flow;
1.308 + }
1.309 +
1.310 +
1.311 +
1.312 +
1.313 + /*
1.314 + Returns the minimum min cut, by a bfs from s in the residual graph.
1.315 + */
1.316 +
1.317 + template<typename CutMap>
1.318 + void mincut(CutMap& M) {
1.319 +
1.320 + std::queue<NodeIt> queue;
1.321 +
1.322 + M.set(s,true);
1.323 + queue.push(s);
1.324 +
1.325 + while (!queue.empty()) {
1.326 + NodeIt w=queue.front();
1.327 + queue.pop();
1.328 +
1.329 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
1.330 + NodeIt v=G.head(e);
1.331 + if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
1.332 + queue.push(v);
1.333 + M.set(v, true);
1.334 + }
1.335 + }
1.336 +
1.337 + for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
1.338 + NodeIt v=G.tail(e);
1.339 + if (!M.get(v) && flow.get(e) > 0 ) {
1.340 + queue.push(v);
1.341 + M.set(v, true);
1.342 + }
1.343 + }
1.344 +
1.345 + }
1.346 +
1.347 + }
1.348 +
1.349 +
1.350 +
1.351 + /*
1.352 + Returns the maximum min cut, by a reverse bfs
1.353 + from t in the residual graph.
1.354 + */
1.355 +
1.356 + template<typename CutMap>
1.357 + void max_mincut(CutMap& M) {
1.358 +
1.359 + std::queue<NodeIt> queue;
1.360 +
1.361 + M.set(t,true);
1.362 + queue.push(t);
1.363 +
1.364 + while (!queue.empty()) {
1.365 + NodeIt w=queue.front();
1.366 + queue.pop();
1.367 +
1.368 + for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
1.369 + NodeIt v=G.tail(e);
1.370 + if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
1.371 + queue.push(v);
1.372 + M.set(v, true);
1.373 + }
1.374 + }
1.375 +
1.376 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
1.377 + NodeIt v=G.head(e);
1.378 + if (!M.get(v) && flow.get(e) > 0 ) {
1.379 + queue.push(v);
1.380 + M.set(v, true);
1.381 + }
1.382 + }
1.383 + }
1.384 +
1.385 + for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) {
1.386 + M.set(v, !M.get(v));
1.387 + }
1.388 +
1.389 + }
1.390 +
1.391 +
1.392 +
1.393 + template<typename CutMap>
1.394 + void min_mincut(CutMap& M) {
1.395 + mincut(M);
1.396 + }
1.397 +
1.398 +
1.399 +
1.400 + };
1.401 +}//namespace marci
1.402 +#endif
1.403 +
1.404 +
1.405 +
1.406 +