1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/jacint/preflow_res.h Fri Apr 23 21:26:32 2004 +0000
1.3 @@ -0,0 +1,497 @@
1.4 +// -*- C++ -*-
1.5 +//The same as preflow.h, using ResGraphWrapper
1.6 +#ifndef HUGO_PREFLOW_RES_H
1.7 +#define HUGO_PREFLOW_RES_H
1.8 +
1.9 +#define H0 20
1.10 +#define H1 1
1.11 +
1.12 +#include <vector>
1.13 +#include <queue>
1.14 +#include <graph_wrapper.h>
1.15 +
1.16 +#include<iostream>
1.17 +
1.18 +namespace hugo {
1.19 +
1.20 + template <typename Graph, typename T,
1.21 + typename CapMap=typename Graph::EdgeMap<T>,
1.22 + typename FlowMap=typename Graph::EdgeMap<T> >
1.23 + class PreflowRes {
1.24 +
1.25 + typedef typename Graph::Node Node;
1.26 + typedef typename Graph::Edge Edge;
1.27 + typedef typename Graph::NodeIt NodeIt;
1.28 + typedef typename Graph::OutEdgeIt OutEdgeIt;
1.29 + typedef typename Graph::InEdgeIt InEdgeIt;
1.30 +
1.31 + const Graph& G;
1.32 + Node s;
1.33 + Node t;
1.34 + const CapMap& capacity;
1.35 + FlowMap& flow;
1.36 + T value;
1.37 + bool constzero;
1.38 +
1.39 + typedef ResGraphWrapper<const Graph, T, CapMap, FlowMap> ResGW;
1.40 + typedef typename ResGW::OutEdgeIt ResOutEdgeIt;
1.41 + typedef typename ResGW::InEdgeIt ResInEdgeIt;
1.42 + typedef typename ResGW::Edge ResEdge;
1.43 +
1.44 + public:
1.45 + PreflowRes(Graph& _G, Node _s, Node _t, CapMap& _capacity,
1.46 + FlowMap& _flow, bool _constzero ) :
1.47 + G(_G), s(_s), t(_t), capacity(_capacity), flow(_flow), constzero(_constzero) {}
1.48 +
1.49 +
1.50 + void run() {
1.51 +
1.52 + ResGW res_graph(G, capacity, flow);
1.53 +
1.54 + value=0; //for the subsequent runs
1.55 +
1.56 + bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd
1.57 + int n=G.nodeNum();
1.58 + int heur0=(int)(H0*n); //time while running 'bound decrease'
1.59 + int heur1=(int)(H1*n); //time while running 'highest label'
1.60 + int heur=heur1; //starting time interval (#of relabels)
1.61 + bool what_heur=1;
1.62 + /*
1.63 + what_heur is 0 in case 'bound decrease'
1.64 + and 1 in case 'highest label'
1.65 + */
1.66 + bool end=false;
1.67 + /*
1.68 + Needed for 'bound decrease', 'true'
1.69 + means no active nodes are above bound b.
1.70 + */
1.71 + int relabel=0;
1.72 + int k=n-2; //bound on the highest level under n containing a node
1.73 + int b=k; //bound on the highest level under n of an active node
1.74 +
1.75 + typename Graph::NodeMap<int> level(G,n);
1.76 + typename Graph::NodeMap<T> excess(G);
1.77 +
1.78 + std::vector<Node> active(n-1,INVALID);
1.79 + typename Graph::NodeMap<Node> next(G,INVALID);
1.80 + //Stack of the active nodes in level i < n.
1.81 + //We use it in both phases.
1.82 +
1.83 + typename Graph::NodeMap<Node> left(G,INVALID);
1.84 + typename Graph::NodeMap<Node> right(G,INVALID);
1.85 + std::vector<Node> level_list(n,INVALID);
1.86 + /*
1.87 + List of the nodes in level i<n.
1.88 + */
1.89 +
1.90 +
1.91 + /*
1.92 + Reverse_bfs from t in the residual graph,
1.93 + to find the starting level.
1.94 + */
1.95 + level.set(t,0);
1.96 + std::queue<Node> bfs_queue;
1.97 + bfs_queue.push(t);
1.98 +
1.99 + while (!bfs_queue.empty()) {
1.100 +
1.101 + Node v=bfs_queue.front();
1.102 + bfs_queue.pop();
1.103 + int l=level[v]+1;
1.104 +
1.105 + ResInEdgeIt e;
1.106 + for(res_graph.first(e,v); res_graph.valid(e);
1.107 + res_graph.next(e)) {
1.108 + Node w=res_graph.tail(e);
1.109 + if ( level[w] == n && w != s ) {
1.110 + bfs_queue.push(w);
1.111 + Node first=level_list[l];
1.112 + if ( G.valid(first) ) left.set(first,w);
1.113 + right.set(w,first);
1.114 + level_list[l]=w;
1.115 + level.set(w, l);
1.116 + }
1.117 + }
1.118 + }
1.119 +
1.120 +
1.121 + if ( !constzero ) {
1.122 + /*
1.123 + Counting the excess
1.124 + */
1.125 + NodeIt v;
1.126 + for(G.first(v); G.valid(v); G.next(v)) {
1.127 + T exc=0;
1.128 +
1.129 + InEdgeIt e;
1.130 + for(G.first(e,v); G.valid(e); G.next(e)) exc+=flow[e];
1.131 + OutEdgeIt f;
1.132 + for(G.first(f,v); G.valid(f); G.next(f)) exc-=flow[e];
1.133 +
1.134 + excess.set(v,exc);
1.135 +
1.136 + //putting the active nodes into the stack
1.137 + int lev=level[v];
1.138 + if ( exc > 0 && lev < n ) {
1.139 + next.set(v,active[lev]);
1.140 + active[lev]=v;
1.141 + }
1.142 + }
1.143 + }
1.144 +
1.145 +
1.146 +
1.147 + //the starting flow
1.148 + ResOutEdgeIt e;
1.149 + for(res_graph.first(e,s); res_graph.valid(e);
1.150 + res_graph.next(e)) {
1.151 + Node w=res_graph.head(e);
1.152 + if ( level[w] < n ) {
1.153 + if ( excess[w] == 0 && w!=t ) {
1.154 + next.set(w,active[level[w]]);
1.155 + active[level[w]]=w;
1.156 + }
1.157 + T rem=res_graph.resCap(e);
1.158 + excess.set(w, excess[w]+rem);
1.159 + res_graph.augment(e, rem );
1.160 + }
1.161 + }
1.162 +
1.163 +
1.164 + /*
1.165 + End of preprocessing
1.166 + */
1.167 +
1.168 +
1.169 +
1.170 + /*
1.171 + Push/relabel on the highest level active nodes.
1.172 + */
1.173 + while ( true ) {
1.174 +
1.175 + if ( b == 0 ) {
1.176 + if ( phase ) break;
1.177 +
1.178 + if ( !what_heur && !end && k > 0 ) {
1.179 + b=k;
1.180 + end=true;
1.181 + } else {
1.182 + phase=1;
1.183 + level.set(s,0);
1.184 + std::queue<Node> bfs_queue;
1.185 + bfs_queue.push(s);
1.186 +
1.187 + while (!bfs_queue.empty()) {
1.188 +
1.189 + Node v=bfs_queue.front();
1.190 + bfs_queue.pop();
1.191 + int l=level[v]+1;
1.192 +
1.193 + ResInEdgeIt e;
1.194 + for(res_graph.first(e,v);
1.195 + res_graph.valid(e); res_graph.next(e)) {
1.196 + Node u=res_graph.tail(e);
1.197 + if ( level[u] >= n ) {
1.198 + bfs_queue.push(u);
1.199 + level.set(u, l);
1.200 + if ( excess[u] > 0 ) {
1.201 + next.set(u,active[l]);
1.202 + active[l]=u;
1.203 + }
1.204 + }
1.205 + }
1.206 +
1.207 + }
1.208 + b=n-2;
1.209 + }
1.210 +
1.211 + }
1.212 +
1.213 +
1.214 + if ( !G.valid(active[b]) ) --b;
1.215 + else {
1.216 + end=false;
1.217 +
1.218 + Node w=active[b];
1.219 + active[b]=next[w];
1.220 + int lev=level[w];
1.221 + T exc=excess[w];
1.222 + int newlevel=n; //bound on the next level of w
1.223 +
1.224 + ResOutEdgeIt e;
1.225 + for(res_graph.first(e,w); res_graph.valid(e); res_graph.next(e)) {
1.226 +
1.227 + Node v=res_graph.head(e);
1.228 + if( lev > level[v] ) {
1.229 + /*Push is allowed now*/
1.230 +
1.231 + if ( excess[v]==0 && v!=t && v!=s ) {
1.232 + int lev_v=level[v];
1.233 + next.set(v,active[lev_v]);
1.234 + active[lev_v]=v;
1.235 + }
1.236 +
1.237 + T remcap=res_graph.resCap(e);
1.238 +
1.239 + if ( remcap >= exc ) {
1.240 + /*A nonsaturating push.*/
1.241 + res_graph.augment(e, exc);
1.242 + excess.set(v, excess[v]+exc);
1.243 + exc=0;
1.244 + break;
1.245 +
1.246 + } else {
1.247 + /*A saturating push.*/
1.248 +
1.249 + res_graph.augment(e, remcap);
1.250 + excess.set(v, excess[v]+remcap);
1.251 + exc-=remcap;
1.252 + }
1.253 + } else if ( newlevel > level[v] ){
1.254 + newlevel = level[v];
1.255 + }
1.256 +
1.257 + }
1.258 +
1.259 + excess.set(w, exc);
1.260 +
1.261 + /*
1.262 + Relabel
1.263 + */
1.264 +
1.265 +
1.266 + if ( exc > 0 ) {
1.267 + //now 'lev' is the old level of w
1.268 +
1.269 + if ( phase ) {
1.270 + level.set(w,++newlevel);
1.271 + next.set(w,active[newlevel]);
1.272 + active[newlevel]=w;
1.273 + b=newlevel;
1.274 + } else {
1.275 + //unlacing starts
1.276 + Node right_n=right[w];
1.277 + Node left_n=left[w];
1.278 +
1.279 + if ( G.valid(right_n) ) {
1.280 + if ( G.valid(left_n) ) {
1.281 + right.set(left_n, right_n);
1.282 + left.set(right_n, left_n);
1.283 + } else {
1.284 + level_list[lev]=right_n;
1.285 + left.set(right_n, INVALID);
1.286 + }
1.287 + } else {
1.288 + if ( G.valid(left_n) ) {
1.289 + right.set(left_n, INVALID);
1.290 + } else {
1.291 + level_list[lev]=INVALID;
1.292 + }
1.293 + }
1.294 + //unlacing ends
1.295 +
1.296 + if ( !G.valid(level_list[lev]) ) {
1.297 +
1.298 + //gapping starts
1.299 + for (int i=lev; i!=k ; ) {
1.300 + Node v=level_list[++i];
1.301 + while ( G.valid(v) ) {
1.302 + level.set(v,n);
1.303 + v=right[v];
1.304 + }
1.305 + level_list[i]=INVALID;
1.306 + if ( !what_heur ) active[i]=INVALID;
1.307 + }
1.308 +
1.309 + level.set(w,n);
1.310 + b=lev-1;
1.311 + k=b;
1.312 + //gapping ends
1.313 +
1.314 + } else {
1.315 +
1.316 + if ( newlevel == n ) level.set(w,n);
1.317 + else {
1.318 + level.set(w,++newlevel);
1.319 + next.set(w,active[newlevel]);
1.320 + active[newlevel]=w;
1.321 + if ( what_heur ) b=newlevel;
1.322 + if ( k < newlevel ) ++k; //now k=newlevel
1.323 + Node first=level_list[newlevel];
1.324 + if ( G.valid(first) ) left.set(first,w);
1.325 + right.set(w,first);
1.326 + left.set(w,INVALID);
1.327 + level_list[newlevel]=w;
1.328 + }
1.329 + }
1.330 +
1.331 +
1.332 + ++relabel;
1.333 + if ( relabel >= heur ) {
1.334 + relabel=0;
1.335 + if ( what_heur ) {
1.336 + what_heur=0;
1.337 + heur=heur0;
1.338 + end=false;
1.339 + } else {
1.340 + what_heur=1;
1.341 + heur=heur1;
1.342 + b=k;
1.343 + }
1.344 + }
1.345 + } //phase 0
1.346 +
1.347 +
1.348 + } // if ( exc > 0 )
1.349 +
1.350 +
1.351 + } // if stack[b] is nonempty
1.352 +
1.353 + } // while(true)
1.354 +
1.355 +
1.356 + value = excess[t];
1.357 + /*Max flow value.*/
1.358 +
1.359 + } //void run()
1.360 +
1.361 +
1.362 +
1.363 +
1.364 +
1.365 + /*
1.366 + Returns the maximum value of a flow.
1.367 + */
1.368 +
1.369 + T flowValue() {
1.370 + return value;
1.371 + }
1.372 +
1.373 +
1.374 + FlowMap Flow() {
1.375 + return flow;
1.376 + }
1.377 +
1.378 +
1.379 +
1.380 + void Flow(FlowMap& _flow ) {
1.381 + NodeIt v;
1.382 + for(G.first(v) ; G.valid(v); G.next(v))
1.383 + _flow.set(v,flow[v]);
1.384 + }
1.385 +
1.386 +
1.387 +
1.388 + /*
1.389 + Returns the minimum min cut, by a bfs from s in the residual graph.
1.390 + */
1.391 +
1.392 + template<typename _CutMap>
1.393 + void minMinCut(_CutMap& M) {
1.394 +
1.395 + std::queue<Node> queue;
1.396 +
1.397 + M.set(s,true);
1.398 + queue.push(s);
1.399 +
1.400 + while (!queue.empty()) {
1.401 + Node w=queue.front();
1.402 + queue.pop();
1.403 +
1.404 + OutEdgeIt e;
1.405 + for(G.first(e,w) ; G.valid(e); G.next(e)) {
1.406 + Node v=G.head(e);
1.407 + if (!M[v] && flow[e] < capacity[e] ) {
1.408 + queue.push(v);
1.409 + M.set(v, true);
1.410 + }
1.411 + }
1.412 +
1.413 + InEdgeIt f;
1.414 + for(G.first(f,w) ; G.valid(f); G.next(f)) {
1.415 + Node v=G.tail(f);
1.416 + if (!M[v] && flow[f] > 0 ) {
1.417 + queue.push(v);
1.418 + M.set(v, true);
1.419 + }
1.420 + }
1.421 + }
1.422 + }
1.423 +
1.424 +
1.425 +
1.426 + /*
1.427 + Returns the maximum min cut, by a reverse bfs
1.428 + from t in the residual graph.
1.429 + */
1.430 +
1.431 + template<typename _CutMap>
1.432 + void maxMinCut(_CutMap& M) {
1.433 +
1.434 + std::queue<Node> queue;
1.435 +
1.436 + M.set(t,true);
1.437 + queue.push(t);
1.438 +
1.439 + while (!queue.empty()) {
1.440 + Node w=queue.front();
1.441 + queue.pop();
1.442 +
1.443 +
1.444 + InEdgeIt e;
1.445 + for(G.first(e,w) ; G.valid(e); G.next(e)) {
1.446 + Node v=G.tail(e);
1.447 + if (!M[v] && flow[e] < capacity[e] ) {
1.448 + queue.push(v);
1.449 + M.set(v, true);
1.450 + }
1.451 + }
1.452 +
1.453 + OutEdgeIt f;
1.454 + for(G.first(f,w) ; G.valid(f); G.next(f)) {
1.455 + Node v=G.head(f);
1.456 + if (!M[v] && flow[f] > 0 ) {
1.457 + queue.push(v);
1.458 + M.set(v, true);
1.459 + }
1.460 + }
1.461 + }
1.462 +
1.463 + NodeIt v;
1.464 + for(G.first(v) ; G.valid(v); G.next(v)) {
1.465 + M.set(v, !M[v]);
1.466 + }
1.467 +
1.468 + }
1.469 +
1.470 +
1.471 +
1.472 + template<typename CutMap>
1.473 + void minCut(CutMap& M) {
1.474 + minMinCut(M);
1.475 + }
1.476 +
1.477 +
1.478 + void reset_target (Node _t) {t=_t;}
1.479 + void reset_source (Node _s) {s=_s;}
1.480 +
1.481 + template<typename _CapMap>
1.482 + void reset_cap (_CapMap _cap) {capacity=_cap;}
1.483 +
1.484 + template<typename _FlowMap>
1.485 + void reset_cap (_FlowMap _flow, bool _constzero) {
1.486 + flow=_flow;
1.487 + constzero=_constzero;
1.488 + }
1.489 +
1.490 +
1.491 +
1.492 + };
1.493 +
1.494 +} //namespace hugo
1.495 +
1.496 +#endif //PREFLOW_H
1.497 +
1.498 +
1.499 +
1.500 +