Flows with test files. The best is preflow.h
1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/jacint/preflow.cc Sat Feb 21 21:01:22 2004 +0000
1.3 @@ -0,0 +1,66 @@
1.4 +#include <iostream>
1.5 +#include <fstream>
1.6 +
1.7 +#include <list_graph.hh>
1.8 +#include <dimacs.hh>
1.9 +#include <preflow.h>
1.10 +#include <time_measure.h>
1.11 +
1.12 +using namespace hugo;
1.13 +
1.14 +// Use a DIMACS max flow file as stdin.
1.15 +// read_dimacs_demo < dimacs_max_flow_file
1.16 +int main(int, char **) {
1.17 + typedef ListGraph::NodeIt NodeIt;
1.18 + typedef ListGraph::EachEdgeIt EachEdgeIt;
1.19 +
1.20 + ListGraph G;
1.21 + NodeIt s, t;
1.22 + ListGraph::EdgeMap<int> cap(G);
1.23 + readDimacsMaxFlow(std::cin, G, s, t, cap);
1.24 +
1.25 + std::cout << "preflow demo ..." << std::endl;
1.26 +
1.27 + double mintime=1000000;
1.28 +
1.29 + for ( int i=1; i!=11; ++i ) {
1.30 + double pre_time=currTime();
1.31 + preflow<ListGraph, int> max_flow_test(G, s, t, cap);
1.32 + double post_time=currTime();
1.33 + if ( mintime > post_time-pre_time ) mintime = post_time-pre_time;
1.34 + }
1.35 +
1.36 + preflow<ListGraph, int> max_flow_test(G, s, t, cap);
1.37 +
1.38 + ListGraph::NodeMap<bool> cut(G);
1.39 + max_flow_test.minCut(cut);
1.40 + int min_cut_value=0;
1.41 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
1.42 + if (cut.get(G.tail(e)) && !cut.get(G.head(e))) min_cut_value+=cap.get(e);
1.43 + }
1.44 +
1.45 + ListGraph::NodeMap<bool> cut1(G);
1.46 + max_flow_test.minMinCut(cut1);
1.47 + int min_min_cut_value=0;
1.48 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
1.49 + if (cut.get(G.tail(e)) && !cut.get(G.head(e)))
1.50 + min_min_cut_value+=cap.get(e);
1.51 + }
1.52 +
1.53 + ListGraph::NodeMap<bool> cut2(G);
1.54 + max_flow_test.maxMinCut(cut2);
1.55 + int max_min_cut_value=0;
1.56 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
1.57 + if (cut2.get(G.tail(e)) && !cut2.get(G.head(e)))
1.58 + max_min_cut_value+=cap.get(e);
1.59 + }
1.60 +
1.61 + std::cout << "min time of 10 runs: " << mintime << " sec"<< std::endl;
1.62 + std::cout << "flow value: "<< max_flow_test.maxFlow() << std::endl;
1.63 + std::cout << "min cut value: "<< min_cut_value << std::endl;
1.64 + std::cout << "min min cut value: "<< min_min_cut_value << std::endl;
1.65 + std::cout << "max min cut value: "<< max_min_cut_value <<
1.66 + std::endl<< std::endl;
1.67 +
1.68 + return 0;
1.69 +}
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/src/work/jacint/preflow.h Sat Feb 21 21:01:22 2004 +0000
2.3 @@ -0,0 +1,532 @@
2.4 +// -*- C++ -*-
2.5 +/*
2.6 +preflow.h
2.7 +by jacint.
2.8 +Heuristics:
2.9 + 2 phase
2.10 + gap
2.11 + list 'level_list' on the nodes on level i implemented by hand
2.12 + stack 'active' on the active nodes on level i implemented by hand
2.13 + runs heuristic 'highest label' for H1*n relabels
2.14 + runs heuristic 'dound decrease' for H0*n relabels, starts with 'highest label'
2.15 +
2.16 +Parameters H0 and H1 are initialized to 20 and 10.
2.17 +
2.18 +The best preflow I could ever write.
2.19 +
2.20 +The constructor runs the algorithm.
2.21 +
2.22 +Members:
2.23 +
2.24 +T maxFlow() : returns the value of a maximum flow
2.25 +
2.26 +T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
2.27 +
2.28 +FlowMap Flow() : returns the fixed maximum flow x
2.29 +
2.30 +void minMinCut(CutMap& M) : sets M to the characteristic vector of the
2.31 + minimum min cut. M should be a map of bools initialized to false.
2.32 +
2.33 +void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
2.34 + maximum min cut. M should be a map of bools initialized to false.
2.35 +
2.36 +void minCut(CutMap& M) : sets M to the characteristic vector of
2.37 + a min cut. M should be a map of bools initialized to false.
2.38 +
2.39 +*/
2.40 +
2.41 +#ifndef PREFLOW_H
2.42 +#define PREFLOW_H
2.43 +
2.44 +#define H0 20
2.45 +#define H1 1
2.46 +
2.47 +#include <vector>
2.48 +#include <queue>
2.49 +
2.50 +namespace hugo {
2.51 +
2.52 + template <typename Graph, typename T,
2.53 + typename FlowMap=typename Graph::EdgeMap<T>,
2.54 + typename CapMap=typename Graph::EdgeMap<T> >
2.55 + class preflow {
2.56 +
2.57 + typedef typename Graph::NodeIt NodeIt;
2.58 + typedef typename Graph::EdgeIt EdgeIt;
2.59 + typedef typename Graph::EachNodeIt EachNodeIt;
2.60 + typedef typename Graph::OutEdgeIt OutEdgeIt;
2.61 + typedef typename Graph::InEdgeIt InEdgeIt;
2.62 +
2.63 + Graph& G;
2.64 + NodeIt s;
2.65 + NodeIt t;
2.66 + FlowMap flow;
2.67 + CapMap& capacity;
2.68 + T value;
2.69 +
2.70 + public:
2.71 +
2.72 + preflow(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity ) :
2.73 + G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity)
2.74 + {
2.75 +
2.76 + bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd
2.77 + int n=G.nodeNum();
2.78 + int heur0=(int)(H0*n); //time while running 'bound decrease'
2.79 + int heur1=(int)(H1*n); //time while running 'highest label'
2.80 + int heur=heur1; //starting time interval (#of relabels)
2.81 + bool what_heur=1;
2.82 + /*
2.83 + what_heur is 0 in case 'bound decrease'
2.84 + and 1 in case 'highest label'
2.85 + */
2.86 + bool end=false;
2.87 + /*
2.88 + Needed for 'bound decrease', 'true'
2.89 + means no active nodes are above bound b.
2.90 + */
2.91 + int relabel=0;
2.92 + int k=n-2; //bound on the highest level under n containing a node
2.93 + int b=k; //bound on the highest level under n of an active node
2.94 +
2.95 + typename Graph::NodeMap<int> level(G,n);
2.96 + typename Graph::NodeMap<T> excess(G);
2.97 +
2.98 + std::vector<NodeIt> active(n);
2.99 + typename Graph::NodeMap<NodeIt> next(G);
2.100 + //Stack of the active nodes in level i < n.
2.101 + //We use it in both phases.
2.102 +
2.103 + typename Graph::NodeMap<NodeIt> left(G);
2.104 + typename Graph::NodeMap<NodeIt> right(G);
2.105 + std::vector<NodeIt> level_list(n);
2.106 + /*
2.107 + List of the nodes in level i<n.
2.108 + */
2.109 +
2.110 + /*Reverse_bfs from t, to find the starting level.*/
2.111 + level.set(t,0);
2.112 + std::queue<NodeIt> bfs_queue;
2.113 + bfs_queue.push(t);
2.114 +
2.115 + while (!bfs_queue.empty()) {
2.116 +
2.117 + NodeIt v=bfs_queue.front();
2.118 + bfs_queue.pop();
2.119 + int l=level.get(v)+1;
2.120 +
2.121 + for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
2.122 + NodeIt w=G.tail(e);
2.123 + if ( level.get(w) == n && w != s ) {
2.124 + bfs_queue.push(w);
2.125 + NodeIt first=level_list[l];
2.126 + if ( first != 0 ) left.set(first,w);
2.127 + right.set(w,first);
2.128 + level_list[l]=w;
2.129 + level.set(w, l);
2.130 + }
2.131 + }
2.132 + }
2.133 +
2.134 + level.set(s,n);
2.135 +
2.136 +
2.137 + /* Starting flow. It is everywhere 0 at the moment. */
2.138 + for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
2.139 + {
2.140 + T c=capacity.get(e);
2.141 + if ( c == 0 ) continue;
2.142 + NodeIt w=G.head(e);
2.143 + if ( level.get(w) < n ) {
2.144 + if ( excess.get(w) == 0 && w!=t ) {
2.145 + next.set(w,active[level.get(w)]);
2.146 + active[level.get(w)]=w;
2.147 + }
2.148 + flow.set(e, c);
2.149 + excess.set(w, excess.get(w)+c);
2.150 + }
2.151 + }
2.152 +
2.153 + /*
2.154 + End of preprocessing
2.155 + */
2.156 +
2.157 +
2.158 +
2.159 + /*
2.160 + Push/relabel on the highest level active nodes.
2.161 + */
2.162 + while ( true ) {
2.163 +
2.164 + if ( b == 0 ) {
2.165 + if ( phase ) break;
2.166 +
2.167 + if ( !what_heur && !end && k > 0 ) {
2.168 + b=k;
2.169 + end=true;
2.170 + } else {
2.171 + phase=1;
2.172 +
2.173 + level.set(s,0);
2.174 + std::queue<NodeIt> bfs_queue;
2.175 + bfs_queue.push(s);
2.176 +
2.177 + while (!bfs_queue.empty()) {
2.178 +
2.179 + NodeIt v=bfs_queue.front();
2.180 + bfs_queue.pop();
2.181 + int l=level.get(v)+1;
2.182 +
2.183 + for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
2.184 + if ( capacity.get(e) == flow.get(e) ) continue;
2.185 + NodeIt u=G.tail(e);
2.186 + if ( level.get(u) >= n ) {
2.187 + bfs_queue.push(u);
2.188 + level.set(u, l);
2.189 + if ( excess.get(u) > 0 ) {
2.190 + next.set(u,active[l]);
2.191 + active[l]=u;
2.192 + }
2.193 + }
2.194 + }
2.195 +
2.196 + for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
2.197 + if ( 0 == flow.get(e) ) continue;
2.198 + NodeIt u=G.head(e);
2.199 + if ( level.get(u) >= n ) {
2.200 + bfs_queue.push(u);
2.201 + level.set(u, l);
2.202 + if ( excess.get(u) > 0 ) {
2.203 + next.set(u,active[l]);
2.204 + active[l]=u;
2.205 + }
2.206 + }
2.207 + }
2.208 + }
2.209 + b=n-2;
2.210 + }
2.211 +
2.212 + }
2.213 +
2.214 +
2.215 + if ( active[b] == 0 ) --b;
2.216 + else {
2.217 + end=false;
2.218 +
2.219 + NodeIt w=active[b];
2.220 + active[b]=next.get(w);
2.221 + int lev=level.get(w);
2.222 + T exc=excess.get(w);
2.223 + int newlevel=n; //bound on the next level of w
2.224 +
2.225 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
2.226 +
2.227 + if ( flow.get(e) == capacity.get(e) ) continue;
2.228 + NodeIt v=G.head(e);
2.229 + //e=wv
2.230 +
2.231 + if( lev > level.get(v) ) {
2.232 + /*Push is allowed now*/
2.233 +
2.234 + if ( excess.get(v)==0 && v!=t && v!=s ) {
2.235 + int lev_v=level.get(v);
2.236 + next.set(v,active[lev_v]);
2.237 + active[lev_v]=v;
2.238 + }
2.239 +
2.240 + T cap=capacity.get(e);
2.241 + T flo=flow.get(e);
2.242 + T remcap=cap-flo;
2.243 +
2.244 + if ( remcap >= exc ) {
2.245 + /*A nonsaturating push.*/
2.246 +
2.247 + flow.set(e, flo+exc);
2.248 + excess.set(v, excess.get(v)+exc);
2.249 + exc=0;
2.250 + break;
2.251 +
2.252 + } else {
2.253 + /*A saturating push.*/
2.254 +
2.255 + flow.set(e, cap);
2.256 + excess.set(v, excess.get(v)+remcap);
2.257 + exc-=remcap;
2.258 + }
2.259 + } else if ( newlevel > level.get(v) ){
2.260 + newlevel = level.get(v);
2.261 + }
2.262 +
2.263 + } //for out edges wv
2.264 +
2.265 +
2.266 + if ( exc > 0 ) {
2.267 + for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
2.268 +
2.269 + if( flow.get(e) == 0 ) continue;
2.270 + NodeIt v=G.tail(e);
2.271 + //e=vw
2.272 +
2.273 + if( lev > level.get(v) ) {
2.274 + /*Push is allowed now*/
2.275 +
2.276 + if ( excess.get(v)==0 && v!=t && v!=s ) {
2.277 + int lev_v=level.get(v);
2.278 + next.set(v,active[lev_v]);
2.279 + active[lev_v]=v;
2.280 + }
2.281 +
2.282 + T flo=flow.get(e);
2.283 +
2.284 + if ( flo >= exc ) {
2.285 + /*A nonsaturating push.*/
2.286 +
2.287 + flow.set(e, flo-exc);
2.288 + excess.set(v, excess.get(v)+exc);
2.289 + exc=0;
2.290 + break;
2.291 + } else {
2.292 + /*A saturating push.*/
2.293 +
2.294 + excess.set(v, excess.get(v)+flo);
2.295 + exc-=flo;
2.296 + flow.set(e,0);
2.297 + }
2.298 + } else if ( newlevel > level.get(v) ) {
2.299 + newlevel = level.get(v);
2.300 + }
2.301 + } //for in edges vw
2.302 +
2.303 + } // if w still has excess after the out edge for cycle
2.304 +
2.305 + excess.set(w, exc);
2.306 +
2.307 + /*
2.308 + Relabel
2.309 + */
2.310 +
2.311 +
2.312 + if ( exc > 0 ) {
2.313 + //now 'lev' is the old level of w
2.314 +
2.315 + if ( phase ) {
2.316 + level.set(w,++newlevel);
2.317 + next.set(w,active[newlevel]);
2.318 + active[newlevel]=w;
2.319 + b=newlevel;
2.320 + } else {
2.321 + //unlacing starts
2.322 + NodeIt right_n=right.get(w);
2.323 + NodeIt left_n=left.get(w);
2.324 +
2.325 + if ( right_n != 0 ) {
2.326 + if ( left_n != 0 ) {
2.327 + right.set(left_n, right_n);
2.328 + left.set(right_n, left_n);
2.329 + } else {
2.330 + level_list[lev]=right_n;
2.331 + left.set(right_n, 0);
2.332 + }
2.333 + } else {
2.334 + if ( left_n != 0 ) {
2.335 + right.set(left_n, 0);
2.336 + } else {
2.337 + level_list[lev]=0;
2.338 +
2.339 + }
2.340 + }
2.341 + //unlacing ends
2.342 +
2.343 + //gapping starts
2.344 + if ( level_list[lev]==0 ) {
2.345 +
2.346 + for (int i=lev; i!=k ; ) {
2.347 + NodeIt v=level_list[++i];
2.348 + while ( v != 0 ) {
2.349 + level.set(v,n);
2.350 + v=right.get(v);
2.351 + }
2.352 + level_list[i]=0;
2.353 + if ( !what_heur ) active[i]=0;
2.354 + }
2.355 +
2.356 + level.set(w,n);
2.357 + b=lev-1;
2.358 + k=b;
2.359 + //gapping ends
2.360 + } else {
2.361 +
2.362 + if ( newlevel == n ) level.set(w,n);
2.363 + else {
2.364 + level.set(w,++newlevel);
2.365 + next.set(w,active[newlevel]);
2.366 + active[newlevel]=w;
2.367 + if ( what_heur ) b=newlevel;
2.368 + if ( k < newlevel ) ++k;
2.369 + NodeIt first=level_list[newlevel];
2.370 + if ( first != 0 ) left.set(first,w);
2.371 + right.set(w,first);
2.372 + left.set(w,0);
2.373 + level_list[newlevel]=w;
2.374 + }
2.375 + }
2.376 +
2.377 +
2.378 + ++relabel;
2.379 + if ( relabel >= heur ) {
2.380 + relabel=0;
2.381 + if ( what_heur ) {
2.382 + what_heur=0;
2.383 + heur=heur0;
2.384 + end=false;
2.385 + } else {
2.386 + what_heur=1;
2.387 + heur=heur1;
2.388 + b=k;
2.389 + }
2.390 + }
2.391 + } //phase 0
2.392 +
2.393 +
2.394 + } // if ( exc > 0 )
2.395 +
2.396 +
2.397 + } // if stack[b] is nonempty
2.398 +
2.399 + } // while(true)
2.400 +
2.401 +
2.402 + value = excess.get(t);
2.403 + /*Max flow value.*/
2.404 +
2.405 + } //void run()
2.406 +
2.407 +
2.408 +
2.409 +
2.410 +
2.411 + /*
2.412 + Returns the maximum value of a flow.
2.413 + */
2.414 +
2.415 + T maxFlow() {
2.416 + return value;
2.417 + }
2.418 +
2.419 +
2.420 +
2.421 + /*
2.422 + For the maximum flow x found by the algorithm,
2.423 + it returns the flow value on edge e, i.e. x(e).
2.424 + */
2.425 +
2.426 + T flowOnEdge(EdgeIt e) {
2.427 + return flow.get(e);
2.428 + }
2.429 +
2.430 +
2.431 +
2.432 + FlowMap Flow() {
2.433 + return flow;
2.434 + }
2.435 +
2.436 +
2.437 +
2.438 + void Flow(FlowMap& _flow ) {
2.439 + for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v)
2.440 + _flow.set(v,flow.get(v));
2.441 + }
2.442 +
2.443 +
2.444 +
2.445 + /*
2.446 + Returns the minimum min cut, by a bfs from s in the residual graph.
2.447 + */
2.448 +
2.449 + template<typename _CutMap>
2.450 + void minMinCut(_CutMap& M) {
2.451 +
2.452 + std::queue<NodeIt> queue;
2.453 +
2.454 + M.set(s,true);
2.455 + queue.push(s);
2.456 +
2.457 + while (!queue.empty()) {
2.458 + NodeIt w=queue.front();
2.459 + queue.pop();
2.460 +
2.461 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
2.462 + NodeIt v=G.head(e);
2.463 + if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
2.464 + queue.push(v);
2.465 + M.set(v, true);
2.466 + }
2.467 + }
2.468 +
2.469 + for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
2.470 + NodeIt v=G.tail(e);
2.471 + if (!M.get(v) && flow.get(e) > 0 ) {
2.472 + queue.push(v);
2.473 + M.set(v, true);
2.474 + }
2.475 + }
2.476 + }
2.477 + }
2.478 +
2.479 +
2.480 +
2.481 + /*
2.482 + Returns the maximum min cut, by a reverse bfs
2.483 + from t in the residual graph.
2.484 + */
2.485 +
2.486 + template<typename _CutMap>
2.487 + void maxMinCut(_CutMap& M) {
2.488 +
2.489 + std::queue<NodeIt> queue;
2.490 +
2.491 + M.set(t,true);
2.492 + queue.push(t);
2.493 +
2.494 + while (!queue.empty()) {
2.495 + NodeIt w=queue.front();
2.496 + queue.pop();
2.497 +
2.498 + for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
2.499 + NodeIt v=G.tail(e);
2.500 + if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
2.501 + queue.push(v);
2.502 + M.set(v, true);
2.503 + }
2.504 + }
2.505 +
2.506 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
2.507 + NodeIt v=G.head(e);
2.508 + if (!M.get(v) && flow.get(e) > 0 ) {
2.509 + queue.push(v);
2.510 + M.set(v, true);
2.511 + }
2.512 + }
2.513 + }
2.514 +
2.515 + for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) {
2.516 + M.set(v, !M.get(v));
2.517 + }
2.518 +
2.519 + }
2.520 +
2.521 +
2.522 +
2.523 + template<typename CutMap>
2.524 + void minCut(CutMap& M) {
2.525 + minMinCut(M);
2.526 + }
2.527 +
2.528 +
2.529 + };
2.530 +}//namespace marci
2.531 +#endif
2.532 +
2.533 +
2.534 +
2.535 +
3.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
3.2 +++ b/src/work/jacint/preflow_hl0.cc Sat Feb 21 21:01:22 2004 +0000
3.3 @@ -0,0 +1,71 @@
3.4 +#include <iostream>
3.5 +#include <fstream>
3.6 +
3.7 +#include <list_graph.hh>
3.8 +#include <dimacs.hh>
3.9 +#include <preflow_hl0.h>
3.10 +#include <time_measure.h>
3.11 +
3.12 +using namespace hugo;
3.13 +
3.14 +// Use a DIMACS max flow file as stdin.
3.15 +// read_dimacs_demo < dimacs_max_flow_file
3.16 +int main(int, char **) {
3.17 + typedef ListGraph::NodeIt NodeIt;
3.18 + typedef ListGraph::EachEdgeIt EachEdgeIt;
3.19 +
3.20 + ListGraph G;
3.21 + NodeIt s, t;
3.22 + ListGraph::EdgeMap<int> cap(G);
3.23 + readDimacsMaxFlow(std::cin, G, s, t, cap);
3.24 +
3.25 + std::cout << "preflow_hl0 demo ..." << std::endl;
3.26 +
3.27 + double mintime=1000000;
3.28 +
3.29 + for ( int i=1; i!=11; ++i ) {
3.30 + double pre_time=currTime();
3.31 + preflow_hl0<ListGraph, int> max_flow_test(G, s, t, cap);
3.32 + double post_time=currTime();
3.33 + if ( mintime > post_time-pre_time ) mintime = post_time-pre_time;
3.34 + }
3.35 +
3.36 + double pre_time=currTime();
3.37 + preflow_hl0<ListGraph, int> max_flow_test(G, s, t, cap);
3.38 + double post_time=currTime();
3.39 +
3.40 + ListGraph::NodeMap<bool> cut(G);
3.41 + max_flow_test.minCut(cut);
3.42 + int min_cut_value=0;
3.43 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
3.44 + if (cut.get(G.tail(e)) && !cut.get(G.head(e))) min_cut_value+=cap.get(e);
3.45 + }
3.46 +
3.47 + ListGraph::NodeMap<bool> cut1(G);
3.48 + max_flow_test.minMinCut(cut1);
3.49 + int min_min_cut_value=0;
3.50 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
3.51 + if (cut.get(G.tail(e)) && !cut.get(G.head(e)))
3.52 + min_min_cut_value+=cap.get(e);
3.53 + }
3.54 +
3.55 + ListGraph::NodeMap<bool> cut2(G);
3.56 + max_flow_test.maxMinCut(cut2);
3.57 + int max_min_cut_value=0;
3.58 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
3.59 + if (cut2.get(G.tail(e)) && !cut2.get(G.head(e)))
3.60 + max_min_cut_value+=cap.get(e);
3.61 + }
3.62 +
3.63 + std::cout << "min time of 10 runs: " << mintime << " sec"<< std::endl;
3.64 + std::cout << "phase 0: " << max_flow_test.time-pre_time
3.65 + << " sec"<< std::endl;
3.66 + std::cout << "phase 1: " << post_time-max_flow_test.time
3.67 + << " sec"<< std::endl;
3.68 + std::cout << "flow value: "<< max_flow_test.maxFlow() << std::endl;
3.69 + std::cout << "min cut value: "<< min_cut_value << std::endl;
3.70 + std::cout << "min min cut value: "<< min_min_cut_value << std::endl;
3.71 + std::cout << "max min cut value: "<< max_min_cut_value << std::endl;
3.72 +
3.73 + return 0;
3.74 +}
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
4.2 +++ b/src/work/jacint/preflow_hl0.h Sat Feb 21 21:01:22 2004 +0000
4.3 @@ -0,0 +1,508 @@
4.4 +// -*- C++ -*-
4.5 +/*
4.6 +preflow_hl0.h
4.7 +by jacint.
4.8 +Heuristics:
4.9 + 2 phase
4.10 + gap
4.11 + list 'level_list' on the nodes on level i implemented by hand
4.12 + stack 'active' on the active nodes on level i implemented by hand
4.13 + bound decrease
4.14 +
4.15 +The bound decrease heuristic behaves unexpectedly well.
4.16 +
4.17 +The constructor runs the algorithm.
4.18 +
4.19 +Members:
4.20 +
4.21 +T maxFlow() : returns the value of a maximum flow
4.22 +
4.23 +T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
4.24 +
4.25 +FlowMap Flow() : returns the fixed maximum flow x
4.26 +
4.27 +void minMinCut(CutMap& M) : sets M to the characteristic vector of the
4.28 + minimum min cut. M should be a map of bools initialized to false.
4.29 +
4.30 +void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
4.31 + maximum min cut. M should be a map of bools initialized to false.
4.32 +
4.33 +
4.34 +void minCut(CutMap& M) : sets M to the characteristic vector of
4.35 + a min cut. M should be a map of bools initialized to false.
4.36 +
4.37 +*/
4.38 +
4.39 +#ifndef PREFLOW_HL0_H
4.40 +#define PREFLOW_HL0_H
4.41 +
4.42 +#include <vector>
4.43 +#include <queue>
4.44 +
4.45 +#include <time_measure.h> //for test
4.46 +
4.47 +namespace hugo {
4.48 +
4.49 + template <typename Graph, typename T,
4.50 + typename FlowMap=typename Graph::EdgeMap<T>,
4.51 + typename CapMap=typename Graph::EdgeMap<T> >
4.52 + class preflow_hl0 {
4.53 +
4.54 + typedef typename Graph::NodeIt NodeIt;
4.55 + typedef typename Graph::EdgeIt EdgeIt;
4.56 + typedef typename Graph::EachNodeIt EachNodeIt;
4.57 + typedef typename Graph::OutEdgeIt OutEdgeIt;
4.58 + typedef typename Graph::InEdgeIt InEdgeIt;
4.59 +
4.60 + Graph& G;
4.61 + NodeIt s;
4.62 + NodeIt t;
4.63 + FlowMap flow;
4.64 + CapMap& capacity;
4.65 + T value;
4.66 +
4.67 + public:
4.68 + double time;
4.69 +
4.70 + preflow_hl0(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity ) :
4.71 + G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) {
4.72 +
4.73 + bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd
4.74 + int n=G.nodeNum();
4.75 + bool end=false;
4.76 + /*
4.77 + 'true' means no active nodes are above bound b.
4.78 + */
4.79 + int k=n-2; //bound on the highest level under n containing a node
4.80 + int b=k; //bound on the highest level under n of an active node
4.81 + /*
4.82 + b is a bound on the highest level of the stack.
4.83 + k is a bound on the highest nonempty level i < n.
4.84 + */
4.85 +
4.86 + typename Graph::NodeMap<int> level(G,n);
4.87 + typename Graph::NodeMap<T> excess(G);
4.88 +
4.89 + std::vector<NodeIt> active(n);
4.90 + typename Graph::NodeMap<NodeIt> next(G);
4.91 + //Stack of the active nodes in level i < n.
4.92 + //We use it in both phases.
4.93 +
4.94 + typename Graph::NodeMap<NodeIt> left(G);
4.95 + typename Graph::NodeMap<NodeIt> right(G);
4.96 + std::vector<NodeIt> level_list(n);
4.97 + /*
4.98 + List of the nodes in level i<n.
4.99 + */
4.100 +
4.101 + /*Reverse_bfs from t, to find the starting level.*/
4.102 + level.set(t,0);
4.103 + std::queue<NodeIt> bfs_queue;
4.104 + bfs_queue.push(t);
4.105 +
4.106 + while (!bfs_queue.empty()) {
4.107 +
4.108 + NodeIt v=bfs_queue.front();
4.109 + bfs_queue.pop();
4.110 + int l=level.get(v)+1;
4.111 +
4.112 + for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
4.113 + NodeIt w=G.tail(e);
4.114 + if ( level.get(w) == n && w != s ) {
4.115 + bfs_queue.push(w);
4.116 + NodeIt first=level_list[l];
4.117 + if ( first != 0 ) left.set(first,w);
4.118 + right.set(w,first);
4.119 + level_list[l]=w;
4.120 + level.set(w, l);
4.121 + }
4.122 + }
4.123 + }
4.124 +
4.125 + level.set(s,n);
4.126 +
4.127 +
4.128 + /* Starting flow. It is everywhere 0 at the moment. */
4.129 + for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
4.130 + {
4.131 + T c=capacity.get(e);
4.132 + if ( c == 0 ) continue;
4.133 + NodeIt w=G.head(e);
4.134 + if ( level.get(w) < n ) {
4.135 + if ( excess.get(w) == 0 && w!=t ) {
4.136 + next.set(w,active[level.get(w)]);
4.137 + active[level.get(w)]=w;
4.138 + }
4.139 + flow.set(e, c);
4.140 + excess.set(w, excess.get(w)+c);
4.141 + }
4.142 + }
4.143 +
4.144 + /*
4.145 + End of preprocessing
4.146 + */
4.147 +
4.148 +
4.149 +
4.150 + /*
4.151 + Push/relabel on the highest level active nodes.
4.152 + */
4.153 + while ( true ) {
4.154 +
4.155 + if ( b == 0 ) {
4.156 + if ( phase ) break;
4.157 +
4.158 + if ( !end && k > 0 ) {
4.159 + b=k;
4.160 + end=true;
4.161 + } else {
4.162 + phase=1;
4.163 + time=currTime();
4.164 + level.set(s,0);
4.165 + std::queue<NodeIt> bfs_queue;
4.166 + bfs_queue.push(s);
4.167 +
4.168 + while (!bfs_queue.empty()) {
4.169 +
4.170 + NodeIt v=bfs_queue.front();
4.171 + bfs_queue.pop();
4.172 + int l=level.get(v)+1;
4.173 +
4.174 + for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
4.175 + if ( capacity.get(e) == flow.get(e) ) continue;
4.176 + NodeIt u=G.tail(e);
4.177 + if ( level.get(u) >= n ) {
4.178 + bfs_queue.push(u);
4.179 + level.set(u, l);
4.180 + if ( excess.get(u) > 0 ) {
4.181 + next.set(u,active[l]);
4.182 + active[l]=u;
4.183 + }
4.184 + }
4.185 + }
4.186 +
4.187 + for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
4.188 + if ( 0 == flow.get(e) ) continue;
4.189 + NodeIt u=G.head(e);
4.190 + if ( level.get(u) >= n ) {
4.191 + bfs_queue.push(u);
4.192 + level.set(u, l);
4.193 + if ( excess.get(u) > 0 ) {
4.194 + next.set(u,active[l]);
4.195 + active[l]=u;
4.196 + }
4.197 + }
4.198 + }
4.199 + }
4.200 + b=n-2;
4.201 + }
4.202 +
4.203 + }
4.204 +
4.205 +
4.206 + if ( active[b] == 0 ) --b;
4.207 + else {
4.208 + end=false;
4.209 +
4.210 + NodeIt w=active[b];
4.211 + active[b]=next.get(w);
4.212 + int lev=level.get(w);
4.213 + T exc=excess.get(w);
4.214 + int newlevel=n; //bound on the next level of w
4.215 +
4.216 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
4.217 +
4.218 + if ( flow.get(e) == capacity.get(e) ) continue;
4.219 + NodeIt v=G.head(e);
4.220 + //e=wv
4.221 +
4.222 + if( lev > level.get(v) ) {
4.223 + /*Push is allowed now*/
4.224 +
4.225 + if ( excess.get(v)==0 && v!=t && v!=s ) {
4.226 + int lev_v=level.get(v);
4.227 + next.set(v,active[lev_v]);
4.228 + active[lev_v]=v;
4.229 + }
4.230 +
4.231 + T cap=capacity.get(e);
4.232 + T flo=flow.get(e);
4.233 + T remcap=cap-flo;
4.234 +
4.235 + if ( remcap >= exc ) {
4.236 + /*A nonsaturating push.*/
4.237 +
4.238 + flow.set(e, flo+exc);
4.239 + excess.set(v, excess.get(v)+exc);
4.240 + exc=0;
4.241 + break;
4.242 +
4.243 + } else {
4.244 + /*A saturating push.*/
4.245 +
4.246 + flow.set(e, cap);
4.247 + excess.set(v, excess.get(v)+remcap);
4.248 + exc-=remcap;
4.249 + }
4.250 + } else if ( newlevel > level.get(v) ){
4.251 + newlevel = level.get(v);
4.252 + }
4.253 +
4.254 + } //for out edges wv
4.255 +
4.256 +
4.257 + if ( exc > 0 ) {
4.258 + for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
4.259 +
4.260 + if( flow.get(e) == 0 ) continue;
4.261 + NodeIt v=G.tail(e);
4.262 + //e=vw
4.263 +
4.264 + if( lev > level.get(v) ) {
4.265 + /*Push is allowed now*/
4.266 +
4.267 + if ( excess.get(v)==0 && v!=t && v!=s ) {
4.268 + int lev_v=level.get(v);
4.269 + next.set(v,active[lev_v]);
4.270 + active[lev_v]=v;
4.271 + }
4.272 +
4.273 + T flo=flow.get(e);
4.274 +
4.275 + if ( flo >= exc ) {
4.276 + /*A nonsaturating push.*/
4.277 +
4.278 + flow.set(e, flo-exc);
4.279 + excess.set(v, excess.get(v)+exc);
4.280 + exc=0;
4.281 + break;
4.282 + } else {
4.283 + /*A saturating push.*/
4.284 +
4.285 + excess.set(v, excess.get(v)+flo);
4.286 + exc-=flo;
4.287 + flow.set(e,0);
4.288 + }
4.289 + } else if ( newlevel > level.get(v) ) {
4.290 + newlevel = level.get(v);
4.291 + }
4.292 + } //for in edges vw
4.293 +
4.294 + } // if w still has excess after the out edge for cycle
4.295 +
4.296 + excess.set(w, exc);
4.297 +
4.298 + /*
4.299 + Relabel
4.300 + */
4.301 +
4.302 +
4.303 + if ( exc > 0 ) {
4.304 + //now 'lev' is the old level of w
4.305 +
4.306 + if ( phase ) {
4.307 + level.set(w,++newlevel);
4.308 + next.set(w,active[newlevel]);
4.309 + active[newlevel]=w;
4.310 + b=newlevel;
4.311 + } else {
4.312 + //unlacing starts
4.313 + NodeIt right_n=right.get(w);
4.314 + NodeIt left_n=left.get(w);
4.315 +
4.316 + if ( right_n != 0 ) {
4.317 + if ( left_n != 0 ) {
4.318 + right.set(left_n, right_n);
4.319 + left.set(right_n, left_n);
4.320 + } else {
4.321 + level_list[lev]=right_n;
4.322 + left.set(right_n, 0);
4.323 + }
4.324 + } else {
4.325 + if ( left_n != 0 ) {
4.326 + right.set(left_n, 0);
4.327 + } else {
4.328 + level_list[lev]=0;
4.329 +
4.330 + }
4.331 + }
4.332 + //unlacing ends
4.333 +
4.334 + //gapping starts
4.335 + if ( level_list[lev]==0 ) {
4.336 +
4.337 + for (int i=lev; i!=k ; ) {
4.338 + NodeIt v=level_list[++i];
4.339 + while ( v != 0 ) {
4.340 + level.set(v,n);
4.341 + v=right.get(v);
4.342 + }
4.343 + level_list[i]=0;
4.344 + active[i]=0;
4.345 + }
4.346 +
4.347 + level.set(w,n);
4.348 + b=lev-1;
4.349 + k=b;
4.350 + //gapping ends
4.351 + } else {
4.352 +
4.353 + if ( newlevel == n ) level.set(w,n);
4.354 + else {
4.355 + level.set(w,++newlevel);
4.356 + next.set(w,active[newlevel]);
4.357 + active[newlevel]=w;
4.358 + if ( k < newlevel ) ++k;
4.359 + NodeIt first=level_list[newlevel];
4.360 + if ( first != 0 ) left.set(first,w);
4.361 + right.set(w,first);
4.362 + left.set(w,0);
4.363 + level_list[newlevel]=w;
4.364 + }
4.365 + }
4.366 +
4.367 +
4.368 + } //phase 0
4.369 +
4.370 +
4.371 + } // if ( exc > 0 )
4.372 +
4.373 +
4.374 + } // if stack[b] is nonempty
4.375 +
4.376 + } // while(true)
4.377 +
4.378 +
4.379 + value = excess.get(t);
4.380 + /*Max flow value.*/
4.381 +
4.382 + } //void run()
4.383 +
4.384 +
4.385 +
4.386 +
4.387 +
4.388 + /*
4.389 + Returns the maximum value of a flow.
4.390 + */
4.391 +
4.392 + T maxFlow() {
4.393 + return value;
4.394 + }
4.395 +
4.396 +
4.397 +
4.398 + /*
4.399 + For the maximum flow x found by the algorithm,
4.400 + it returns the flow value on edge e, i.e. x(e).
4.401 + */
4.402 +
4.403 + T flowOnEdge(EdgeIt e) {
4.404 + return flow.get(e);
4.405 + }
4.406 +
4.407 +
4.408 +
4.409 + FlowMap Flow() {
4.410 + return flow;
4.411 + }
4.412 +
4.413 +
4.414 +
4.415 + void Flow(FlowMap& _flow ) {
4.416 + for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v)
4.417 + _flow.set(v,flow.get(v));
4.418 + }
4.419 +
4.420 +
4.421 +
4.422 + /*
4.423 + Returns the minimum min cut, by a bfs from s in the residual graph.
4.424 + */
4.425 +
4.426 + template<typename _CutMap>
4.427 + void minMinCut(_CutMap& M) {
4.428 +
4.429 + std::queue<NodeIt> queue;
4.430 +
4.431 + M.set(s,true);
4.432 + queue.push(s);
4.433 +
4.434 + while (!queue.empty()) {
4.435 + NodeIt w=queue.front();
4.436 + queue.pop();
4.437 +
4.438 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
4.439 + NodeIt v=G.head(e);
4.440 + if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
4.441 + queue.push(v);
4.442 + M.set(v, true);
4.443 + }
4.444 + }
4.445 +
4.446 + for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
4.447 + NodeIt v=G.tail(e);
4.448 + if (!M.get(v) && flow.get(e) > 0 ) {
4.449 + queue.push(v);
4.450 + M.set(v, true);
4.451 + }
4.452 + }
4.453 + }
4.454 + }
4.455 +
4.456 +
4.457 +
4.458 + /*
4.459 + Returns the maximum min cut, by a reverse bfs
4.460 + from t in the residual graph.
4.461 + */
4.462 +
4.463 + template<typename _CutMap>
4.464 + void maxMinCut(_CutMap& M) {
4.465 +
4.466 + std::queue<NodeIt> queue;
4.467 +
4.468 + M.set(t,true);
4.469 + queue.push(t);
4.470 +
4.471 + while (!queue.empty()) {
4.472 + NodeIt w=queue.front();
4.473 + queue.pop();
4.474 +
4.475 + for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
4.476 + NodeIt v=G.tail(e);
4.477 + if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
4.478 + queue.push(v);
4.479 + M.set(v, true);
4.480 + }
4.481 + }
4.482 +
4.483 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
4.484 + NodeIt v=G.head(e);
4.485 + if (!M.get(v) && flow.get(e) > 0 ) {
4.486 + queue.push(v);
4.487 + M.set(v, true);
4.488 + }
4.489 + }
4.490 + }
4.491 +
4.492 + for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) {
4.493 + M.set(v, !M.get(v));
4.494 + }
4.495 +
4.496 + }
4.497 +
4.498 +
4.499 +
4.500 + template<typename _CutMap>
4.501 + void minCut(_CutMap& M) {
4.502 + minMinCut(M);
4.503 + }
4.504 +
4.505 + };
4.506 +}//namespace marci
4.507 +#endif
4.508 +
4.509 +
4.510 +
4.511 +
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
5.2 +++ b/src/work/jacint/preflow_hl1.cc Sat Feb 21 21:01:22 2004 +0000
5.3 @@ -0,0 +1,73 @@
5.4 +#include <iostream>
5.5 +#include <fstream>
5.6 +
5.7 +#include <list_graph.hh>
5.8 +#include <dimacs.hh>
5.9 +#include <preflow_hl1.h>
5.10 +#include <time_measure.h>
5.11 +
5.12 +using namespace hugo;
5.13 +
5.14 +// Use a DIMACS max flow file as stdin.
5.15 +// read_dimacs_demo < dimacs_max_flow_file
5.16 +int main(int, char **) {
5.17 + typedef ListGraph::NodeIt NodeIt;
5.18 + typedef ListGraph::EachEdgeIt EachEdgeIt;
5.19 +
5.20 + ListGraph G;
5.21 + NodeIt s, t;
5.22 + ListGraph::EdgeMap<int> cap(G);
5.23 + readDimacsMaxFlow(std::cin, G, s, t, cap);
5.24 +
5.25 + std::cout << "preflow_hl1 demo ..." << std::endl;
5.26 +
5.27 + double mintime=1000000;
5.28 +
5.29 + for ( int i=1; i!=11; ++i ) {
5.30 + double pre_time=currTime();
5.31 + preflow_hl1<ListGraph, int> max_flow_test(G, s, t, cap);
5.32 + double post_time=currTime();
5.33 + if ( mintime > post_time-pre_time ) mintime = post_time-pre_time;
5.34 + }
5.35 +
5.36 + double pre_time=currTime();
5.37 + preflow_hl1<ListGraph, int> max_flow_test(G, s, t, cap);
5.38 + double post_time=currTime();
5.39 +
5.40 + ListGraph::NodeMap<bool> cut(G);
5.41 + max_flow_test.minCut(cut);
5.42 + int min_cut_value=0;
5.43 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
5.44 + if (cut.get(G.tail(e)) && !cut.get(G.head(e))) min_cut_value+=cap.get(e);
5.45 + }
5.46 +
5.47 + ListGraph::NodeMap<bool> cut1(G);
5.48 + max_flow_test.minMinCut(cut1);
5.49 + int min_min_cut_value=0;
5.50 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
5.51 + if (cut.get(G.tail(e)) && !cut.get(G.head(e)))
5.52 + min_min_cut_value+=cap.get(e);
5.53 + }
5.54 +
5.55 + ListGraph::NodeMap<bool> cut2(G);
5.56 + max_flow_test.maxMinCut(cut2);
5.57 + int max_min_cut_value=0;
5.58 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
5.59 + if (cut2.get(G.tail(e)) && !cut2.get(G.head(e)))
5.60 + max_min_cut_value+=cap.get(e);
5.61 + }
5.62 +
5.63 + std::cout << "min time of 10 runs: " << mintime << " sec"<< std::endl;
5.64 + std::cout << "phase 0: " << max_flow_test.time-pre_time
5.65 + << " sec"<< std::endl;
5.66 + std::cout << "phase 1: " << post_time-max_flow_test.time
5.67 + << " sec"<< std::endl;
5.68 + std::cout << "flow value: "<< max_flow_test.maxFlow() << std::endl;
5.69 + std::cout << "min cut value: "<< min_cut_value << std::endl;
5.70 + std::cout << "min min cut value: "<< min_min_cut_value << std::endl;
5.71 + std::cout << "max min cut value: "<< max_min_cut_value <<
5.72 + std::endl<< std::endl;
5.73 +
5.74 +
5.75 + return 0;
5.76 +}
6.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
6.2 +++ b/src/work/jacint/preflow_hl2.cc Sat Feb 21 21:01:22 2004 +0000
6.3 @@ -0,0 +1,65 @@
6.4 +#include <iostream>
6.5 +#include <fstream>
6.6 +
6.7 +#include <list_graph.hh>
6.8 +#include <dimacs.hh>
6.9 +#include <preflow_hl2.h>
6.10 +#include <time_measure.h>
6.11 +
6.12 +using namespace hugo;
6.13 +
6.14 +// Use a DIMACS max flow file as stdin.
6.15 +// read_dimacs_demo < dimacs_max_flow_file
6.16 +int main(int, char **) {
6.17 + typedef ListGraph::NodeIt NodeIt;
6.18 + typedef ListGraph::EachEdgeIt EachEdgeIt;
6.19 +
6.20 + ListGraph G;
6.21 + NodeIt s, t;
6.22 + ListGraph::EdgeMap<int> cap(G);
6.23 + readDimacsMaxFlow(std::cin, G, s, t, cap);
6.24 +
6.25 + std::cout << "preflow_hl2 demo ..." << std::endl;
6.26 +
6.27 + double mintime=1000000;
6.28 +
6.29 + for ( int i=1; i!=11; ++i ) {
6.30 + double pre_time=currTime();
6.31 + preflow_hl2<ListGraph, int> max_flow_test(G, s, t, cap);
6.32 + double post_time=currTime();
6.33 + if ( mintime > post_time-pre_time ) mintime = post_time-pre_time;
6.34 + }
6.35 +
6.36 + preflow_hl2<ListGraph, int> max_flow_test(G, s, t, cap);
6.37 +
6.38 + ListGraph::NodeMap<bool> cut(G);
6.39 + max_flow_test.minCut(cut);
6.40 + int min_cut_value=0;
6.41 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
6.42 + if (cut.get(G.tail(e)) && !cut.get(G.head(e))) min_cut_value+=cap.get(e);
6.43 + }
6.44 +
6.45 + ListGraph::NodeMap<bool> cut1(G);
6.46 + max_flow_test.minMinCut(cut1);
6.47 + int min_min_cut_value=0;
6.48 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
6.49 + if (cut.get(G.tail(e)) && !cut.get(G.head(e)))
6.50 + min_min_cut_value+=cap.get(e);
6.51 + }
6.52 +
6.53 + ListGraph::NodeMap<bool> cut2(G);
6.54 + max_flow_test.maxMinCut(cut2);
6.55 + int max_min_cut_value=0;
6.56 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
6.57 + if (cut2.get(G.tail(e)) && !cut2.get(G.head(e)))
6.58 + max_min_cut_value+=cap.get(e);
6.59 + }
6.60 +
6.61 + std::cout << "min time of 10 runs: " << mintime << " sec"<< std::endl;
6.62 + std::cout << "flow value: "<< max_flow_test.maxFlow() << std::endl;
6.63 + std::cout << "min cut value: "<< min_cut_value << std::endl;
6.64 + std::cout << "min min cut value: "<< min_min_cut_value << std::endl;
6.65 + std::cout << "max min cut value: "<< max_min_cut_value << std::endl;
6.66 +
6.67 + return 0;
6.68 +}
7.1 --- a/src/work/jacint/preflow_hl2.h Fri Feb 20 22:01:02 2004 +0000
7.2 +++ b/src/work/jacint/preflow_hl2.h Sat Feb 21 21:01:22 2004 +0000
7.3 @@ -3,31 +3,35 @@
7.4 preflow_hl2.h
7.5 by jacint.
7.6 Runs the highest label variant of the preflow push algorithm with
7.7 -running time O(n^2\sqrt(m)), with the 'empty level' and with the
7.8 -heuristic that the bound b on the active nodes is not increased
7.9 -only when b=0, when we put b=2*n-2.
7.10 +running time O(n^2\sqrt(m)).
7.11
7.12 -'A' is a parameter for the empty_level heuristic
7.13 +Heuristics:
7.14
7.15 -Member functions:
7.16 + gap: we iterate through the nodes for finding the nodes above
7.17 + the gap and under level n. So it is quite slow.
7.18 + numb: we maintain the number of nodes in level i.
7.19 + highest label
7.20 +
7.21 +'A' is a parameter for the gap, we only upgrade the nodes to level n,
7.22 + if the gap is under A*n.
7.23
7.24 -void run() : runs the algorithm
7.25 +The constructor runs the algorithm.
7.26
7.27 - The following functions should be used after run() was already run.
7.28 +Members:
7.29
7.30 -T maxflow() : returns the value of a maximum flow
7.31 +T maxFlow() : returns the value of a maximum flow
7.32
7.33 -T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
7.34 +T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
7.35
7.36 -FlowMap allflow() : returns the fixed maximum flow x
7.37 +FlowMap Flow() : returns the fixed maximum flow x
7.38
7.39 -void mincut(CutMap& M) : sets M to the characteristic vector of a
7.40 +void minCut(CutMap& M) : sets M to the characteristic vector of a
7.41 minimum cut. M should be a map of bools initialized to false.
7.42
7.43 -void min_mincut(CutMap& M) : sets M to the characteristic vector of the
7.44 +void minMinCut(CutMap& M) : sets M to the characteristic vector of the
7.45 minimum min cut. M should be a map of bools initialized to false.
7.46
7.47 -void max_mincut(CutMap& M) : sets M to the characteristic vector of the
7.48 +void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
7.49 maximum min cut. M should be a map of bools initialized to false.
7.50
7.51 */
7.52 @@ -35,7 +39,7 @@
7.53 #ifndef PREFLOW_HL2_H
7.54 #define PREFLOW_HL2_H
7.55
7.56 -#define A 1
7.57 +#define A .9
7.58
7.59 #include <vector>
7.60 #include <stack>
7.61 @@ -44,8 +48,8 @@
7.62 namespace hugo {
7.63
7.64 template <typename Graph, typename T,
7.65 - typename FlowMap=typename Graph::EdgeMap<T>, typename CapMap=typename Graph::EdgeMap<T>,
7.66 - typename IntMap=typename Graph::NodeMap<int>, typename TMap=typename Graph::NodeMap<T> >
7.67 + typename FlowMap=typename Graph::EdgeMap<T>,
7.68 + typename CapMap=typename Graph::EdgeMap<T> >
7.69 class preflow_hl2 {
7.70
7.71 typedef typename Graph::NodeIt NodeIt;
7.72 @@ -64,22 +68,17 @@
7.73 public:
7.74
7.75 preflow_hl2(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
7.76 - G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) { }
7.77 + G(_G), s(_s), t(_t), flow(_G), capacity(_capacity) {
7.78
7.79 -
7.80 - void run() {
7.81 -
7.82 - bool no_end=true;
7.83 int n=G.nodeNum();
7.84 int b=n-2;
7.85 /*
7.86 b is a bound on the highest level of an active node.
7.87 - In the beginning it is at most n-2.
7.88 */
7.89
7.90 - IntMap level(G,n);
7.91 - TMap excess(G);
7.92 -
7.93 + typename Graph::NodeMap<int> level(G,n);
7.94 + typename Graph::NodeMap<T> excess(G);
7.95 +
7.96 std::vector<int> numb(n);
7.97 /*
7.98 The number of nodes on level i < n. It is
7.99 @@ -110,7 +109,7 @@
7.100 }
7.101 }
7.102 }
7.103 -
7.104 +
7.105 level.set(s,n);
7.106
7.107
7.108 @@ -127,38 +126,28 @@
7.109 excess.set(w, excess.get(w)+c);
7.110 }
7.111 }
7.112 -
7.113 +
7.114 /*
7.115 End of preprocessing
7.116 */
7.117
7.118
7.119 -
7.120 /*
7.121 Push/relabel on the highest level active nodes.
7.122 - */
7.123 + */
7.124 /*While there exists an active node.*/
7.125 while (b) {
7.126 - if ( stack[b].empty() ) {
7.127 - if ( b==1 ) {
7.128 - if ( !no_end ) break;
7.129 - else {
7.130 - b=2*n-2;
7.131 - no_end=false;
7.132 - }
7.133 - }
7.134 + if ( stack[b].empty() ) {
7.135 --b;
7.136 - } else {
7.137 -
7.138 - no_end=true;
7.139 -
7.140 - NodeIt w=stack[b].top(); //w is a highest label active node.
7.141 - stack[b].pop();
7.142 - int lev=level.get(w);
7.143 - int exc=excess.get(w);
7.144 - int newlevel=2*n; //In newlevel we bound the next level of w.
7.145 -
7.146 - // if ( level.get(w) < n ) { //Nem tudom ez mukodik-e
7.147 + continue;
7.148 + }
7.149 +
7.150 + NodeIt w=stack[b].top();
7.151 + stack[b].pop();
7.152 + int lev=level.get(w);
7.153 + T exc=excess.get(w);
7.154 + int newlevel=2*n; //In newlevel we bound the next level of w.
7.155 +
7.156 for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
7.157
7.158 if ( flow.get(e) == capacity.get(e) ) continue;
7.159 @@ -172,9 +161,9 @@
7.160 stack[level.get(v)].push(v);
7.161 /*v becomes active.*/
7.162
7.163 - int cap=capacity.get(e);
7.164 - int flo=flow.get(e);
7.165 - int remcap=cap-flo;
7.166 + T cap=capacity.get(e);
7.167 + T flo=flow.get(e);
7.168 + T remcap=cap-flo;
7.169
7.170 if ( remcap >= exc ) {
7.171 /*A nonsaturating push.*/
7.172 @@ -210,7 +199,7 @@
7.173 stack[level.get(v)].push(v);
7.174 /*v becomes active.*/
7.175
7.176 - int flo=flow.get(e);
7.177 + T flo=flow.get(e);
7.178
7.179 if ( flo >= exc ) {
7.180 /*A nonsaturating push.*/
7.181 @@ -231,42 +220,43 @@
7.182 } //for in edges vw
7.183
7.184 } // if w still has excess after the out edge for cycle
7.185 -
7.186 - excess.set(w, exc);
7.187 +
7.188 + excess.set(w, exc);
7.189 +
7.190 +
7.191 +
7.192 +
7.193 + /*
7.194 + Relabel
7.195 + */
7.196 +
7.197 + if ( exc > 0 ) {
7.198 + //now 'lev' is the old level of w
7.199 + level.set(w,++newlevel);
7.200
7.201 -
7.202 - /*
7.203 - Relabel
7.204 - */
7.205 + if ( lev < n ) {
7.206 + --numb[lev];
7.207 +
7.208 + if ( !numb[lev] && lev < A*n ) { //If the level of w gets empty.
7.209 +
7.210 + for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
7.211 + if (level.get(v) > lev && level.get(v) < n ) level.set(v,n);
7.212 + }
7.213 + for (int i=lev+1 ; i!=n ; ++i) numb[i]=0;
7.214 + if ( newlevel < n ) newlevel=n;
7.215 + } else {
7.216 + if ( newlevel < n ) ++numb[newlevel];
7.217 + }
7.218 + }
7.219
7.220 - if ( exc > 0 ) {
7.221 - //now 'lev' is the old level of w
7.222 - level.set(w,++newlevel);
7.223 -
7.224 - if ( lev < n ) {
7.225 - --numb[lev];
7.226 -
7.227 - if ( !numb[lev] && lev < A*n ) { //If the level of w gets empty.
7.228 -
7.229 - for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
7.230 - if (level.get(v) > lev && level.get(v) < n ) level.set(v,n);
7.231 - }
7.232 - for (int i=lev+1 ; i!=n ; ++i) numb[i]=0;
7.233 - if ( newlevel < n ) newlevel=n;
7.234 - } else {
7.235 - if ( newlevel < n ) ++numb[newlevel];
7.236 - }
7.237 - }
7.238 -
7.239 - stack[newlevel].push(w);
7.240 -
7.241 - }
7.242 -
7.243 - } // if stack[b] is nonempty
7.244 -
7.245 + stack[newlevel].push(w);
7.246 + b=newlevel;
7.247 +
7.248 + }
7.249 +
7.250 } // while(b)
7.251 -
7.252 -
7.253 +
7.254 +
7.255 value = excess.get(t);
7.256 /*Max flow value.*/
7.257
7.258 @@ -281,17 +271,18 @@
7.259 Returns the maximum value of a flow.
7.260 */
7.261
7.262 - T maxflow() {
7.263 + T maxFlow() {
7.264 return value;
7.265 }
7.266
7.267
7.268
7.269 /*
7.270 - For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e).
7.271 + For the maximum flow x found by the algorithm,
7.272 + it returns the flow value on edge e, i.e. x(e).
7.273 */
7.274
7.275 - T flowonedge(EdgeIt e) {
7.276 + T flowOnEdge(const EdgeIt e) {
7.277 return flow.get(e);
7.278 }
7.279
7.280 @@ -301,7 +292,7 @@
7.281 Returns the maximum flow x found by the algorithm.
7.282 */
7.283
7.284 - FlowMap allflow() {
7.285 + FlowMap Flow() {
7.286 return flow;
7.287 }
7.288
7.289 @@ -313,7 +304,7 @@
7.290 */
7.291
7.292 template<typename CutMap>
7.293 - void mincut(CutMap& M) {
7.294 + void minCut(CutMap& M) {
7.295
7.296 std::queue<NodeIt> queue;
7.297
7.298 @@ -323,7 +314,7 @@
7.299 while (!queue.empty()) {
7.300 NodeIt w=queue.front();
7.301 queue.pop();
7.302 -
7.303 +
7.304 for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
7.305 NodeIt v=G.head(e);
7.306 if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
7.307 @@ -338,10 +329,9 @@
7.308 queue.push(v);
7.309 M.set(v, true);
7.310 }
7.311 - }
7.312 + }
7.313
7.314 }
7.315 -
7.316 }
7.317
7.318
7.319 @@ -352,7 +342,7 @@
7.320 */
7.321
7.322 template<typename CutMap>
7.323 - void max_mincut(CutMap& M) {
7.324 + void maxMinCut(CutMap& M) {
7.325
7.326 std::queue<NodeIt> queue;
7.327
7.328 @@ -389,14 +379,14 @@
7.329
7.330
7.331 template<typename CutMap>
7.332 - void min_mincut(CutMap& M) {
7.333 - mincut(M);
7.334 + void minMinCut(CutMap& M) {
7.335 + minCut(M);
7.336 }
7.337
7.338
7.339
7.340 };
7.341 -}//namespace hugo
7.342 +}//namespace marci
7.343 #endif
7.344
7.345
8.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
8.2 +++ b/src/work/jacint/preflow_hl3.cc Sat Feb 21 21:01:22 2004 +0000
8.3 @@ -0,0 +1,72 @@
8.4 +#include <iostream>
8.5 +#include <fstream>
8.6 +
8.7 +#include <list_graph.hh>
8.8 +#include <dimacs.hh>
8.9 +#include <preflow_hl3.h>
8.10 +#include <time_measure.h>
8.11 +
8.12 +using namespace hugo;
8.13 +
8.14 +// Use a DIMACS max flow file as stdin.
8.15 +// read_dimacs_demo < dimacs_max_flow_file
8.16 +int main(int, char **) {
8.17 + typedef ListGraph::NodeIt NodeIt;
8.18 + typedef ListGraph::EachEdgeIt EachEdgeIt;
8.19 +
8.20 + ListGraph G;
8.21 + NodeIt s, t;
8.22 + ListGraph::EdgeMap<int> cap(G);
8.23 + readDimacsMaxFlow(std::cin, G, s, t, cap);
8.24 +
8.25 + std::cout << "preflow_hl3 demo ..." << std::endl;
8.26 +
8.27 + double mintime=1000000;
8.28 +
8.29 + for ( int i=1; i!=11; ++i ) {
8.30 + double pre_time=currTime();
8.31 + preflow_hl3<ListGraph, int> max_flow_test(G, s, t, cap);
8.32 + double post_time=currTime();
8.33 + if ( mintime > post_time-pre_time ) mintime = post_time-pre_time;
8.34 + }
8.35 +
8.36 + double pre_time=currTime();
8.37 + preflow_hl3<ListGraph, int> max_flow_test(G, s, t, cap);
8.38 + double post_time=currTime();
8.39 +
8.40 + ListGraph::NodeMap<bool> cut(G);
8.41 + max_flow_test.minCut(cut);
8.42 + int min_cut_value=0;
8.43 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
8.44 + if (cut.get(G.tail(e)) && !cut.get(G.head(e))) min_cut_value+=cap.get(e);
8.45 + }
8.46 +
8.47 + ListGraph::NodeMap<bool> cut1(G);
8.48 + max_flow_test.minMinCut(cut1);
8.49 + int min_min_cut_value=0;
8.50 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
8.51 + if (cut.get(G.tail(e)) && !cut.get(G.head(e)))
8.52 + min_min_cut_value+=cap.get(e);
8.53 + }
8.54 +
8.55 + ListGraph::NodeMap<bool> cut2(G);
8.56 + max_flow_test.maxMinCut(cut2);
8.57 + int max_min_cut_value=0;
8.58 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
8.59 + if (cut2.get(G.tail(e)) && !cut2.get(G.head(e)))
8.60 + max_min_cut_value+=cap.get(e);
8.61 + }
8.62 +
8.63 + std::cout << "min time of 10 runs: " << mintime << " sec"<< std::endl;
8.64 + std::cout << "phase 0: " << max_flow_test.time-pre_time
8.65 + << " sec"<< std::endl;
8.66 + std::cout << "phase 1: " << post_time-max_flow_test.time
8.67 + << " sec"<< std::endl;
8.68 + std::cout << "flow value: "<< max_flow_test.maxFlow() << std::endl;
8.69 + std::cout << "min cut value: "<< min_cut_value << std::endl;
8.70 + std::cout << "min min cut value: "<< min_min_cut_value << std::endl;
8.71 + std::cout << "max min cut value: "<< max_min_cut_value <<
8.72 + std::endl<< std::endl;
8.73 +
8.74 + return 0;
8.75 +}
9.1 --- a/src/work/jacint/preflow_hl3.h Fri Feb 20 22:01:02 2004 +0000
9.2 +++ b/src/work/jacint/preflow_hl3.h Sat Feb 21 21:01:22 2004 +0000
9.3 @@ -2,35 +2,31 @@
9.4 /*
9.5 preflow_hl3.h
9.6 by jacint.
9.7 -Runs the highest label variant of the preflow push algorithm with
9.8 -running time O(n^2\sqrt(m)), with the felszippantos 'empty level'
9.9 -and with the two-phase heuristic: if there is no active node of
9.10 -level at most n, then we go into phase 1, do a bfs
9.11 -from s, and flow the excess back to s.
9.12
9.13 -In phase 1 we shift everything downwards by n.
9.14 +Heuristics:
9.15
9.16 -'A' is a parameter for the empty_level heuristic
9.17 + suck gap : if there is a gap, then we put all
9.18 + nodes reachable from the relabelled node to level n
9.19 + 2 phase
9.20 + highest label
9.21
9.22 -Member functions:
9.23 +The constructor runs the algorithm.
9.24
9.25 -void run() : runs the algorithm
9.26 +Members:
9.27
9.28 - The following functions should be used after run() was already run.
9.29 +T maxFlow() : returns the value of a maximum flow
9.30
9.31 -T maxflow() : returns the value of a maximum flow
9.32 +T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
9.33
9.34 -T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
9.35 +FlowMap Flow() : returns the fixed maximum flow x
9.36
9.37 -FlowMap allflow() : returns the fixed maximum flow x
9.38 -
9.39 -void mincut(CutMap& M) : sets M to the characteristic vector of a
9.40 +void minCut(CutMap& M) : sets M to the characteristic vector of a
9.41 minimum cut. M should be a map of bools initialized to false.
9.42
9.43 -void min_mincut(CutMap& M) : sets M to the characteristic vector of the
9.44 +void minMinCut(CutMap& M) : sets M to the characteristic vector of the
9.45 minimum min cut. M should be a map of bools initialized to false.
9.46
9.47 -void max_mincut(CutMap& M) : sets M to the characteristic vector of the
9.48 +void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
9.49 maximum min cut. M should be a map of bools initialized to false.
9.50
9.51 */
9.52 @@ -38,17 +34,17 @@
9.53 #ifndef PREFLOW_HL3_H
9.54 #define PREFLOW_HL3_H
9.55
9.56 -#define A 1
9.57 -
9.58 #include <vector>
9.59 #include <stack>
9.60 #include <queue>
9.61
9.62 +#include <time_measure.h> //for test
9.63 +
9.64 namespace hugo {
9.65
9.66 template <typename Graph, typename T,
9.67 - typename FlowMap=typename Graph::EdgeMap<T>, typename CapMap=typename Graph::EdgeMap<T>,
9.68 - typename IntMap=typename Graph::NodeMap<int>, typename TMap=typename Graph::NodeMap<T> >
9.69 + typename FlowMap=typename Graph::EdgeMap<T>,
9.70 + typename CapMap=typename Graph::EdgeMap<T> >
9.71 class preflow_hl3 {
9.72
9.73 typedef typename Graph::NodeIt NodeIt;
9.74 @@ -66,12 +62,11 @@
9.75
9.76 public:
9.77
9.78 + double time; //for test
9.79 +
9.80 preflow_hl3(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
9.81 - G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) { }
9.82 -
9.83 -
9.84 - void run() {
9.85 -
9.86 + G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) {
9.87 +
9.88 bool phase=0;
9.89 int n=G.nodeNum();
9.90 int b=n-2;
9.91 @@ -80,8 +75,8 @@
9.92 In the beginning it is at most n-2.
9.93 */
9.94
9.95 - IntMap level(G,n);
9.96 - TMap excess(G);
9.97 + typename Graph::NodeMap<int> level(G,n);
9.98 + typename Graph::NodeMap<T> excess(G);
9.99
9.100 std::vector<int> numb(n);
9.101 /*
9.102 @@ -148,7 +143,8 @@
9.103 if ( b == 0 ) {
9.104 if ( phase ) break;
9.105 phase=1;
9.106 -
9.107 + time=currTime();
9.108 +
9.109 level.set(s,0);
9.110
9.111 std::queue<NodeIt> bfs_queue;
9.112 @@ -187,11 +183,11 @@
9.113 if ( stack[b].empty() ) --b;
9.114 else {
9.115
9.116 - NodeIt w=stack[b].top(); //w is a highest label active node.
9.117 + NodeIt w=stack[b].top();
9.118 stack[b].pop();
9.119 int lev=level.get(w);
9.120 - int exc=excess.get(w);
9.121 - int newlevel=n; //In newlevel we bound the next level of w.
9.122 + T exc=excess.get(w);
9.123 + int newlevel=n; //bound on the next level of w.
9.124
9.125 for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
9.126
9.127 @@ -206,9 +202,9 @@
9.128 stack[level.get(v)].push(v);
9.129 /*v becomes active.*/
9.130
9.131 - int cap=capacity.get(e);
9.132 - int flo=flow.get(e);
9.133 - int remcap=cap-flo;
9.134 + T cap=capacity.get(e);
9.135 + T flo=flow.get(e);
9.136 + T remcap=cap-flo;
9.137
9.138 if ( remcap >= exc ) {
9.139 /*A nonsaturating push.*/
9.140 @@ -244,7 +240,7 @@
9.141 stack[level.get(v)].push(v);
9.142 /*v becomes active.*/
9.143
9.144 - int flo=flow.get(e);
9.145 + T flo=flow.get(e);
9.146
9.147 if ( flo >= exc ) {
9.148 /*A nonsaturating push.*/
9.149 @@ -349,17 +345,18 @@
9.150 Returns the maximum value of a flow.
9.151 */
9.152
9.153 - T maxflow() {
9.154 + T maxFlow() {
9.155 return value;
9.156 }
9.157
9.158
9.159
9.160 /*
9.161 - For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e).
9.162 + For the maximum flow x found by the algorithm,
9.163 + it returns the flow value on edge e, i.e. x(e).
9.164 */
9.165
9.166 - T flowonedge(EdgeIt e) {
9.167 + T flowOnEdge(EdgeIt e) {
9.168 return flow.get(e);
9.169 }
9.170
9.171 @@ -369,7 +366,7 @@
9.172 Returns the maximum flow x found by the algorithm.
9.173 */
9.174
9.175 - FlowMap allflow() {
9.176 + FlowMap Flow() {
9.177 return flow;
9.178 }
9.179
9.180 @@ -381,7 +378,7 @@
9.181 */
9.182
9.183 template<typename CutMap>
9.184 - void mincut(CutMap& M) {
9.185 + void minCut(CutMap& M) {
9.186
9.187 std::queue<NodeIt> queue;
9.188
9.189 @@ -420,7 +417,7 @@
9.190 */
9.191
9.192 template<typename CutMap>
9.193 - void max_mincut(CutMap& M) {
9.194 + void maxMinCut(CutMap& M) {
9.195
9.196 std::queue<NodeIt> queue;
9.197
9.198 @@ -457,14 +454,14 @@
9.199
9.200
9.201 template<typename CutMap>
9.202 - void min_mincut(CutMap& M) {
9.203 - mincut(M);
9.204 + void minMinCut(CutMap& M) {
9.205 + minCut(M);
9.206 }
9.207
9.208
9.209
9.210 };
9.211 -}//namespace hugo
9.212 +}//namespace
9.213 #endif
9.214
9.215
10.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
10.2 +++ b/src/work/jacint/preflow_hl4.cc Sat Feb 21 21:01:22 2004 +0000
10.3 @@ -0,0 +1,76 @@
10.4 +#include <iostream>
10.5 +#include <fstream>
10.6 +
10.7 +#include <list_graph.hh>
10.8 +#include <dimacs.hh>
10.9 +#include <preflow_hl4.h>
10.10 +#include <time_measure.h>
10.11 +
10.12 +using namespace hugo;
10.13 +
10.14 +//Note, that preflow_hl4.h has a member NodeMap<bool> cut,
10.15 +//the construction of which slowing the running time by 1-10%.
10.16 +
10.17 +// Use a DIMACS max flow file as stdin.
10.18 +// read_dimacs_demo < dimacs_max_flow_file
10.19 +int main(int, char **) {
10.20 + typedef ListGraph::NodeIt NodeIt;
10.21 + typedef ListGraph::EachEdgeIt EachEdgeIt;
10.22 +
10.23 + ListGraph G;
10.24 + NodeIt s, t;
10.25 + ListGraph::EdgeMap<int> cap(G);
10.26 + readDimacsMaxFlow(std::cin, G, s, t, cap);
10.27 +
10.28 + std::cout << "preflow_hl4 demo ..." << std::endl;
10.29 +
10.30 + double mintime=1000000;
10.31 +
10.32 + for ( int i=1; i!=11; ++i ) {
10.33 + double pre_time=currTime();
10.34 + preflow_hl4<ListGraph, int> max_flow_test(G, s, t, cap);
10.35 + double post_time=currTime();
10.36 + if ( mintime > post_time-pre_time ) mintime = post_time-pre_time;
10.37 + }
10.38 +
10.39 + double pre_time=currTime();
10.40 + preflow_hl4<ListGraph, int> max_flow_test(G, s, t, cap);
10.41 + double post_time=currTime();
10.42 +
10.43 + ListGraph::NodeMap<bool> cut(G);
10.44 + max_flow_test.minCut(cut);
10.45 + int min_cut_value=0;
10.46 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
10.47 + if (cut.get(G.tail(e)) && !cut.get(G.head(e))) min_cut_value+=cap.get(e);
10.48 + }
10.49 +
10.50 + ListGraph::NodeMap<bool> cut1(G);
10.51 + max_flow_test.minMinCut(cut1);
10.52 + int min_min_cut_value=0;
10.53 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
10.54 + if (cut.get(G.tail(e)) && !cut.get(G.head(e)))
10.55 + min_min_cut_value+=cap.get(e);
10.56 + }
10.57 +
10.58 + ListGraph::NodeMap<bool> cut2(G);
10.59 + max_flow_test.maxMinCut(cut2);
10.60 + int max_min_cut_value=0;
10.61 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
10.62 + if (cut2.get(G.tail(e)) && !cut2.get(G.head(e)))
10.63 + max_min_cut_value+=cap.get(e);
10.64 + }
10.65 +
10.66 + std::cout << "min time of 10 runs: " << mintime << " sec"<< std::endl;
10.67 + std::cout << "phase 0: " << max_flow_test.time-pre_time
10.68 + << " sec"<< std::endl;
10.69 + std::cout << "phase 1: " << post_time-max_flow_test.time
10.70 + << " sec"<< std::endl;
10.71 + std::cout << "flow value: "<< max_flow_test.maxFlow() << std::endl;
10.72 + std::cout << "min cut value: "<< min_cut_value << std::endl;
10.73 + std::cout << "min min cut value: "<< min_min_cut_value << std::endl;
10.74 + std::cout << "max min cut value: "<< max_min_cut_value <<
10.75 + std::endl<< std::endl;
10.76 +
10.77 +
10.78 + return 0;
10.79 +}
11.1 --- a/src/work/jacint/preflow_hl4.h Fri Feb 20 22:01:02 2004 +0000
11.2 +++ b/src/work/jacint/preflow_hl4.h Sat Feb 21 21:01:22 2004 +0000
11.3 @@ -1,53 +1,67 @@
11.4 // -*- C++ -*-
11.5 /*
11.6 -preflow_hl4.h
11.7 +preflow_h5.h
11.8 by jacint.
11.9 -Runs the two phase highest label preflow push algorithm. In phase 0
11.10 -we maintain in a list the nodes in level i < n, and we maintain a
11.11 -bound k on the max level i < n containing a node, so we can do
11.12 -the gap heuristic fast. Phase 1 is the same. (The algorithm is the
11.13 -same as preflow.hl3, the only diff is that here we use the gap
11.14 -heuristic with the list of the nodes on level i, and not a bfs form the
11.15 -upgraded node.)
11.16 +Heuristics:
11.17 + 2 phase
11.18 + gap
11.19 + list 'level_list' on the nodes on level i implemented by hand
11.20 + highest label
11.21 + relevel: in phase 0, after BFS*n relabels, it runs a reverse
11.22 + bfs from t in the res graph to relevel the nodes reachable from t.
11.23 + BFS is initialized to 20
11.24
11.25 -In phase 1 we shift everything downwards by n.
11.26 +Due to the last heuristic, this algorithm is quite fast on very
11.27 +sparse graphs, but relatively bad on even the dense graphs.
11.28
11.29 -Member functions:
11.30 +'NodeMap<bool> cut' is a member, in this way we can count it fast, after
11.31 +the algorithm was run.
11.32
11.33 -void run() : runs the algorithm
11.34 +The constructor runs the algorithm.
11.35
11.36 - The following functions should be used after run() was already run.
11.37 +Members:
11.38
11.39 -T maxflow() : returns the value of a maximum flow
11.40 +T maxFlow() : returns the value of a maximum flow
11.41
11.42 -T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
11.43 +T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
11.44
11.45 -FlowMap allflow() : returns a maximum flow
11.46 +FlowMap Flow() : returns the fixed maximum flow x
11.47
11.48 -void allflow(FlowMap& _flow ) : returns a maximum flow
11.49 +void Flow(FlowMap& _flow ) : returns the fixed maximum flow x
11.50
11.51 -void mincut(CutMap& M) : sets M to the characteristic vector of a
11.52 - minimum cut. M should be a map of bools initialized to false.
11.53 -
11.54 -void min_mincut(CutMap& M) : sets M to the characteristic vector of the
11.55 +void minMinCut(CutMap& M) : sets M to the characteristic vector of the
11.56 minimum min cut. M should be a map of bools initialized to false.
11.57
11.58 -void max_mincut(CutMap& M) : sets M to the characteristic vector of the
11.59 +void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
11.60 maximum min cut. M should be a map of bools initialized to false.
11.61
11.62 +void minCut(CutMap& M) : fast function, sets M to the characteristic
11.63 + vector of a minimum cut.
11.64 +
11.65 +Different member from the other preflow_hl-s (here we have a member
11.66 +'NodeMap<bool> cut').
11.67 +
11.68 +CutMap minCut() : fast function, giving the characteristic
11.69 + vector of a minimum cut.
11.70 +
11.71 +
11.72 */
11.73
11.74 #ifndef PREFLOW_HL4_H
11.75 #define PREFLOW_HL4_H
11.76
11.77 +#define BFS 20
11.78 +
11.79 #include <vector>
11.80 -#include <stack>
11.81 #include <queue>
11.82
11.83 +#include <time_measure.h> //for test
11.84 +
11.85 namespace hugo {
11.86
11.87 template <typename Graph, typename T,
11.88 typename FlowMap=typename Graph::EdgeMap<T>,
11.89 + typename CutMap=typename Graph::NodeMap<bool>,
11.90 typename CapMap=typename Graph::EdgeMap<T> >
11.91 class preflow_hl4 {
11.92
11.93 @@ -62,18 +76,21 @@
11.94 NodeIt t;
11.95 FlowMap flow;
11.96 CapMap& capacity;
11.97 + CutMap cut;
11.98 T value;
11.99
11.100 public:
11.101
11.102 + double time;
11.103 +
11.104 preflow_hl4(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
11.105 - G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) { }
11.106 + G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity),
11.107 + cut(G, false) {
11.108
11.109 -
11.110 - void run() {
11.111 -
11.112 - bool phase=0;
11.113 + bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd
11.114 int n=G.nodeNum();
11.115 + int relabel=0;
11.116 + int heur=(int)BFS*n;
11.117 int k=n-2;
11.118 int b=k;
11.119 /*
11.120 @@ -83,7 +100,9 @@
11.121
11.122 typename Graph::NodeMap<int> level(G,n);
11.123 typename Graph::NodeMap<T> excess(G);
11.124 - std::vector<std::stack<NodeIt> > stack(n);
11.125 +
11.126 + std::vector<NodeIt> active(n);
11.127 + typename Graph::NodeMap<NodeIt> next(G);
11.128 //Stack of the active nodes in level i < n.
11.129 //We use it in both phases.
11.130
11.131 @@ -107,7 +126,7 @@
11.132
11.133 for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
11.134 NodeIt w=G.tail(e);
11.135 - if ( level.get(w) == n ) {
11.136 + if ( level.get(w) == n && w !=s ) {
11.137 bfs_queue.push(w);
11.138 NodeIt first=level_list[l];
11.139 if ( first != 0 ) left.set(first,w);
11.140 @@ -128,7 +147,10 @@
11.141 if ( c == 0 ) continue;
11.142 NodeIt w=G.head(e);
11.143 if ( level.get(w) < n ) {
11.144 - if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w);
11.145 + if ( excess.get(w) == 0 && w!=t ) {
11.146 + next.set(w,active[level.get(w)]);
11.147 + active[level.get(w)]=w;
11.148 + }
11.149 flow.set(e, c);
11.150 excess.set(w, excess.get(w)+c);
11.151 }
11.152 @@ -151,6 +173,13 @@
11.153 the residual graph.
11.154 */
11.155 phase=1;
11.156 +
11.157 + //Now have a min cut.
11.158 + for( EachNodeIt v=G.template first<EachNodeIt>();
11.159 + v.valid(); ++v)
11.160 + if (level.get(v) >= n ) cut.set(v,true);
11.161 +
11.162 + time=currTime();
11.163 level.set(s,0);
11.164 std::queue<NodeIt> bfs_queue;
11.165 bfs_queue.push(s);
11.166 @@ -167,7 +196,10 @@
11.167 if ( level.get(u) >= n ) {
11.168 bfs_queue.push(u);
11.169 level.set(u, l);
11.170 - if ( excess.get(u) > 0 ) stack[l].push(u);
11.171 + if ( excess.get(u) > 0 ) {
11.172 + next.set(u,active[l]);
11.173 + active[l]=u;
11.174 + }
11.175 }
11.176 }
11.177
11.178 @@ -177,7 +209,10 @@
11.179 if ( level.get(u) >= n ) {
11.180 bfs_queue.push(u);
11.181 level.set(u, l);
11.182 - if ( excess.get(u) > 0 ) stack[l].push(u);
11.183 + if ( excess.get(u) > 0 ) {
11.184 + next.set(u,active[l]);
11.185 + active[l]=u;
11.186 + }
11.187 }
11.188 }
11.189 }
11.190 @@ -185,14 +220,14 @@
11.191 }
11.192
11.193
11.194 - if ( stack[b].empty() ) --b;
11.195 + if ( active[b] == 0 ) --b;
11.196 else {
11.197
11.198 - NodeIt w=stack[b].top(); //w is a highest label active node.
11.199 - stack[b].pop();
11.200 + NodeIt w=active[b];
11.201 + active[b]=next.get(w);
11.202 int lev=level.get(w);
11.203 T exc=excess.get(w);
11.204 - int newlevel=n; //In newlevel we bound the next level of w.
11.205 + int newlevel=n; //bound on the next level of w.
11.206
11.207 for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
11.208
11.209 @@ -203,9 +238,11 @@
11.210 if( lev > level.get(v) ) {
11.211 /*Push is allowed now*/
11.212
11.213 - if ( excess.get(v)==0 && v!=t && v!=s )
11.214 - stack[level.get(v)].push(v);
11.215 - /*v becomes active.*/
11.216 + if ( excess.get(v)==0 && v!=t && v!=s ) {
11.217 + int lev_v=level.get(v);
11.218 + next.set(v,active[lev_v]);
11.219 + active[lev_v]=v;
11.220 + }
11.221
11.222 T cap=capacity.get(e);
11.223 T flo=flow.get(e);
11.224 @@ -243,9 +280,11 @@
11.225 if( lev > level.get(v) ) {
11.226 /*Push is allowed now*/
11.227
11.228 - if ( excess.get(v)==0 && v!=t && v!=s )
11.229 - stack[level.get(v)].push(v);
11.230 - /*v becomes active.*/
11.231 + if ( excess.get(v)==0 && v!=t && v!=s ) {
11.232 + int lev_v=level.get(v);
11.233 + next.set(v,active[lev_v]);
11.234 + active[lev_v]=v;
11.235 + }
11.236
11.237 T flo=flow.get(e);
11.238
11.239 @@ -281,7 +320,8 @@
11.240
11.241 if ( phase ) {
11.242 level.set(w,++newlevel);
11.243 - stack[newlevel].push(w);
11.244 + next.set(w,active[newlevel]);
11.245 + active[newlevel]=w;
11.246 b=newlevel;
11.247 } else {
11.248 //unlacing
11.249 @@ -303,6 +343,7 @@
11.250 level_list[lev]=0;
11.251 }
11.252 }
11.253 + //unlacing ends
11.254
11.255
11.256 if ( level_list[lev]==0 ) {
11.257 @@ -328,7 +369,8 @@
11.258 } else {
11.259
11.260 level.set(w,++newlevel);
11.261 - stack[newlevel].push(w);
11.262 + next.set(w,active[newlevel]);
11.263 + active[newlevel]=w;
11.264 b=newlevel;
11.265 if ( k < newlevel ) ++k;
11.266 NodeIt first=level_list[newlevel];
11.267 @@ -338,9 +380,78 @@
11.268 level_list[newlevel]=w;
11.269 }
11.270 }
11.271 +
11.272 + ++relabel;
11.273 + if ( relabel >= heur ) {
11.274 + relabel=0;
11.275 + b=n-2;
11.276 + k=b;
11.277 +
11.278 + for ( int i=1; i!=n; ++i ) {
11.279 + active[i]=0;
11.280 + level_list[i]=0;
11.281 + }
11.282 +
11.283 + //bfs from t in the res graph to relevel the nodes
11.284 + for( EachNodeIt v=G.template first<EachNodeIt>();
11.285 + v.valid(); ++v) level.set(v,n);
11.286 +
11.287 + level.set(t,0);
11.288 + std::queue<NodeIt> bfs_queue;
11.289 + bfs_queue.push(t);
11.290 +
11.291 + while (!bfs_queue.empty()) {
11.292 +
11.293 + NodeIt v=bfs_queue.front();
11.294 + bfs_queue.pop();
11.295 + int l=level.get(v)+1;
11.296 +
11.297 + for(InEdgeIt e=G.template first<InEdgeIt>(v);
11.298 + e.valid(); ++e) {
11.299 + if ( capacity.get(e) == flow.get(e) ) continue;
11.300 + NodeIt u=G.tail(e);
11.301 + if ( level.get(u) == n ) {
11.302 + bfs_queue.push(u);
11.303 + level.set(u, l);
11.304 + if ( excess.get(u) > 0 ) {
11.305 + next.set(u,active[l]);
11.306 + active[l]=u;
11.307 + }
11.308 + NodeIt first=level_list[l];
11.309 + if ( first != 0 ) left.set(first,w);
11.310 + right.set(w,first);
11.311 + left.set(w,0);
11.312 + level_list[l]=w;
11.313 + }
11.314 + }
11.315 +
11.316 +
11.317 + for(OutEdgeIt e=G.template first<OutEdgeIt>(v);
11.318 + e.valid(); ++e) {
11.319 + if ( 0 == flow.get(e) ) continue;
11.320 + NodeIt u=G.head(e);
11.321 + if ( level.get(u) == n ) {
11.322 + bfs_queue.push(u);
11.323 + level.set(u, l);
11.324 + if ( excess.get(u) > 0 ) {
11.325 + next.set(u,active[l]);
11.326 + active[l]=u;
11.327 + }
11.328 + NodeIt first=level_list[l];
11.329 + if ( first != 0 ) left.set(first,w);
11.330 + right.set(w,first);
11.331 + left.set(w,0);
11.332 + level_list[l]=w;
11.333 + }
11.334 + }
11.335 + }
11.336 +
11.337 + level.set(s,n);
11.338 + }
11.339 +
11.340 } //phase 0
11.341 } // if ( exc > 0 )
11.342 -
11.343 +
11.344
11.345 } // if stack[b] is nonempty
11.346
11.347 @@ -361,7 +472,7 @@
11.348 Returns the maximum value of a flow.
11.349 */
11.350
11.351 - T maxflow() {
11.352 + T maxFlow() {
11.353 return value;
11.354 }
11.355
11.356 @@ -371,19 +482,19 @@
11.357 For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e).
11.358 */
11.359
11.360 - T flowonedge(EdgeIt e) {
11.361 + T flowOnEdge(EdgeIt e) {
11.362 return flow.get(e);
11.363 }
11.364
11.365
11.366
11.367 - FlowMap allflow() {
11.368 + FlowMap Flow() {
11.369 return flow;
11.370 }
11.371
11.372
11.373
11.374 - void allflow(FlowMap& _flow ) {
11.375 + void Flow(FlowMap& _flow ) {
11.376 for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v)
11.377 _flow.set(v,flow.get(v));
11.378 }
11.379 @@ -394,8 +505,8 @@
11.380 Returns the minimum min cut, by a bfs from s in the residual graph.
11.381 */
11.382
11.383 - template<typename CutMap>
11.384 - void mincut(CutMap& M) {
11.385 + template<typename _CutMap>
11.386 + void minMinCut(_CutMap& M) {
11.387
11.388 std::queue<NodeIt> queue;
11.389
11.390 @@ -433,8 +544,8 @@
11.391 from t in the residual graph.
11.392 */
11.393
11.394 - template<typename CutMap>
11.395 - void max_mincut(CutMap& M) {
11.396 + template<typename _CutMap>
11.397 + void maxMinCut(_CutMap& M) {
11.398
11.399 std::queue<NodeIt> queue;
11.400
11.401 @@ -469,16 +580,21 @@
11.402 }
11.403
11.404
11.405 -
11.406 - template<typename CutMap>
11.407 - void min_mincut(CutMap& M) {
11.408 - mincut(M);
11.409 + template<typename _CutMap>
11.410 + void minCut(_CutMap& M) {
11.411 + for( EachNodeIt v=G.template first<EachNodeIt>();
11.412 + v.valid(); ++v)
11.413 + M.set(v, cut.get(v));
11.414 }
11.415
11.416 +
11.417 + CutMap minCut() {
11.418 + return cut;
11.419 + }
11.420
11.421
11.422 };
11.423 -}//namespace hugo
11.424 +}//namespace marci
11.425 #endif
11.426
11.427
12.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
12.2 +++ b/src/work/jacint/preflow_max_flow.cc Sat Feb 21 21:01:22 2004 +0000
12.3 @@ -0,0 +1,40 @@
12.4 +#include <iostream>
12.5 +#include <fstream>
12.6 +
12.7 +#include <list_graph.hh>
12.8 +#include <dimacs.hh>
12.9 +#include <preflow_max_flow.h>
12.10 +#include <time_measure.h>
12.11 +
12.12 +using namespace hugo;
12.13 +
12.14 +// Use a DIMACS max flow file as stdin.
12.15 +// read_dimacs_demo < dimacs_max_flow_file
12.16 +int main(int, char **) {
12.17 + typedef ListGraph::NodeIt NodeIt;
12.18 + typedef ListGraph::EachEdgeIt EachEdgeIt;
12.19 +typedef ListGraph::EachNodeIt EachNodeIt;
12.20 +
12.21 + ListGraph G;
12.22 + NodeIt s, t;
12.23 + ListGraph::EdgeMap<int> cap(G);
12.24 + readDimacsMaxFlow(std::cin, G, s, t, cap);
12.25 +
12.26 + std::cout << "preflow_max_flow demo ..." << std::endl;
12.27 +
12.28 + double pre_time=currTime();
12.29 + preflow_max_flow<ListGraph, int> max_flow_test(G, s, t, cap);
12.30 + ListGraph::NodeMap<bool> cut=max_flow_test.minCut();
12.31 + double post_time=currTime();
12.32 +
12.33 + int cut_value=0;
12.34 + for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
12.35 + if (cut.get(G.tail(e)) && !cut.get(G.head(e))) cut_value+=cap.get(e);
12.36 + }
12.37 +
12.38 + std::cout << "elapsed time: " << post_time-pre_time << " sec"<< std::endl;
12.39 + std::cout << "flow value: "<< max_flow_test.maxFlow() << std::endl;
12.40 + std::cout << "cut value: "<< cut_value << std::endl;
12.41 +
12.42 + return 0;
12.43 +}
13.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
13.2 +++ b/src/work/jacint/preflow_max_flow.h Sat Feb 21 21:01:22 2004 +0000
13.3 @@ -0,0 +1,358 @@
13.4 +// -*- C++ -*-
13.5 +/*
13.6 +preflow_max_flow.h
13.7 +by jacint.
13.8 +Runs the first phase of preflow.h
13.9 +
13.10 +The constructor runs the algorithm.
13.11 +
13.12 +Members:
13.13 +
13.14 +T maxFlow() : returns the value of a maximum flow
13.15 +
13.16 +CutMap minCut() : returns the characteristic vector of a min cut.
13.17 +*/
13.18 +
13.19 +#ifndef PREFLOW_MAX_FLOW_H
13.20 +#define PREFLOW_MAX_FLOW_H
13.21 +
13.22 +#define H0 20
13.23 +#define H1 1
13.24 +
13.25 +#include <vector>
13.26 +#include <queue>
13.27 +
13.28 +namespace hugo {
13.29 +
13.30 + template <typename Graph, typename T,
13.31 + typename FlowMap=typename Graph::EdgeMap<T>,
13.32 + typename CapMap=typename Graph::EdgeMap<T>,
13.33 + typename CutMap=typename Graph::NodeMap<bool> >
13.34 + class preflow_max_flow {
13.35 +
13.36 + typedef typename Graph::NodeIt NodeIt;
13.37 + typedef typename Graph::EdgeIt EdgeIt;
13.38 + typedef typename Graph::EachNodeIt EachNodeIt;
13.39 + typedef typename Graph::OutEdgeIt OutEdgeIt;
13.40 + typedef typename Graph::InEdgeIt InEdgeIt;
13.41 +
13.42 + Graph& G;
13.43 + NodeIt s;
13.44 + NodeIt t;
13.45 + FlowMap flow;
13.46 + CapMap& capacity;
13.47 + CutMap cut;
13.48 + T value;
13.49 +
13.50 + public:
13.51 +
13.52 + preflow_max_flow(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity ) :
13.53 + G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity), cut(_G, false)
13.54 + {
13.55 +
13.56 + int n=G.nodeNum();
13.57 + int heur0=(int)(H0*n); //time while running 'bound decrease'
13.58 + int heur1=(int)(H1*n); //time while running 'highest label'
13.59 + int heur=heur1; //starting time interval (#of relabels)
13.60 + bool what_heur=1;
13.61 + /*
13.62 + what_heur is 0 in case 'bound decrease'
13.63 + and 1 in case 'highest label'
13.64 + */
13.65 + bool end=false;
13.66 + /*
13.67 + Needed for 'bound decrease', 'true'
13.68 + means no active nodes are above bound b.
13.69 + */
13.70 + int relabel=0;
13.71 + int k=n-2; //bound on the highest level under n containing a node
13.72 + int b=k; //bound on the highest level under n of an active node
13.73 +
13.74 + typename Graph::NodeMap<int> level(G,n);
13.75 + typename Graph::NodeMap<T> excess(G);
13.76 +
13.77 + std::vector<NodeIt> active(n);
13.78 + typename Graph::NodeMap<NodeIt> next(G);
13.79 + //Stack of the active nodes in level i < n.
13.80 + //We use it in both phases.
13.81 +
13.82 + typename Graph::NodeMap<NodeIt> left(G);
13.83 + typename Graph::NodeMap<NodeIt> right(G);
13.84 + std::vector<NodeIt> level_list(n);
13.85 + /*
13.86 + List of the nodes in level i<n.
13.87 + */
13.88 +
13.89 + /*Reverse_bfs from t, to find the starting level.*/
13.90 + level.set(t,0);
13.91 + std::queue<NodeIt> bfs_queue;
13.92 + bfs_queue.push(t);
13.93 +
13.94 + while (!bfs_queue.empty()) {
13.95 +
13.96 + NodeIt v=bfs_queue.front();
13.97 + bfs_queue.pop();
13.98 + int l=level.get(v)+1;
13.99 +
13.100 + for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
13.101 + NodeIt w=G.tail(e);
13.102 + if ( level.get(w) == n && w != s ) {
13.103 + bfs_queue.push(w);
13.104 + NodeIt first=level_list[l];
13.105 + if ( first != 0 ) left.set(first,w);
13.106 + right.set(w,first);
13.107 + level_list[l]=w;
13.108 + level.set(w, l);
13.109 + }
13.110 + }
13.111 + }
13.112 +
13.113 + level.set(s,n);
13.114 +
13.115 +
13.116 + /* Starting flow. It is everywhere 0 at the moment. */
13.117 + for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
13.118 + {
13.119 + T c=capacity.get(e);
13.120 + if ( c == 0 ) continue;
13.121 + NodeIt w=G.head(e);
13.122 + if ( level.get(w) < n ) {
13.123 + if ( excess.get(w) == 0 && w!=t ) {
13.124 + next.set(w,active[level.get(w)]);
13.125 + active[level.get(w)]=w;
13.126 + }
13.127 + flow.set(e, c);
13.128 + excess.set(w, excess.get(w)+c);
13.129 + }
13.130 + }
13.131 +
13.132 + /*
13.133 + End of preprocessing
13.134 + */
13.135 +
13.136 +
13.137 +
13.138 + /*
13.139 + Push/relabel on the highest level active nodes.
13.140 + */
13.141 + while ( true ) {
13.142 +
13.143 + if ( b == 0 ) {
13.144 + if ( !what_heur && !end && k > 0 ) {
13.145 + b=k;
13.146 + end=true;
13.147 + } else break;
13.148 + }
13.149 +
13.150 +
13.151 + if ( active[b] == 0 ) --b;
13.152 + else {
13.153 + end=false;
13.154 +
13.155 + NodeIt w=active[b];
13.156 + active[b]=next.get(w);
13.157 + int lev=level.get(w);
13.158 + T exc=excess.get(w);
13.159 + int newlevel=n; //bound on the next level of w
13.160 +
13.161 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
13.162 +
13.163 + if ( flow.get(e) == capacity.get(e) ) continue;
13.164 + NodeIt v=G.head(e);
13.165 + //e=wv
13.166 +
13.167 + if( lev > level.get(v) ) {
13.168 + /*Push is allowed now*/
13.169 +
13.170 + if ( excess.get(v)==0 && v!=t && v!=s ) {
13.171 + int lev_v=level.get(v);
13.172 + next.set(v,active[lev_v]);
13.173 + active[lev_v]=v;
13.174 + }
13.175 +
13.176 + T cap=capacity.get(e);
13.177 + T flo=flow.get(e);
13.178 + T remcap=cap-flo;
13.179 +
13.180 + if ( remcap >= exc ) {
13.181 + /*A nonsaturating push.*/
13.182 +
13.183 + flow.set(e, flo+exc);
13.184 + excess.set(v, excess.get(v)+exc);
13.185 + exc=0;
13.186 + break;
13.187 +
13.188 + } else {
13.189 + /*A saturating push.*/
13.190 +
13.191 + flow.set(e, cap);
13.192 + excess.set(v, excess.get(v)+remcap);
13.193 + exc-=remcap;
13.194 + }
13.195 + } else if ( newlevel > level.get(v) ){
13.196 + newlevel = level.get(v);
13.197 + }
13.198 +
13.199 + } //for out edges wv
13.200 +
13.201 +
13.202 + if ( exc > 0 ) {
13.203 + for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
13.204 +
13.205 + if( flow.get(e) == 0 ) continue;
13.206 + NodeIt v=G.tail(e);
13.207 + //e=vw
13.208 +
13.209 + if( lev > level.get(v) ) {
13.210 + /*Push is allowed now*/
13.211 +
13.212 + if ( excess.get(v)==0 && v!=t && v!=s ) {
13.213 + int lev_v=level.get(v);
13.214 + next.set(v,active[lev_v]);
13.215 + active[lev_v]=v;
13.216 + }
13.217 +
13.218 + T flo=flow.get(e);
13.219 +
13.220 + if ( flo >= exc ) {
13.221 + /*A nonsaturating push.*/
13.222 +
13.223 + flow.set(e, flo-exc);
13.224 + excess.set(v, excess.get(v)+exc);
13.225 + exc=0;
13.226 + break;
13.227 + } else {
13.228 + /*A saturating push.*/
13.229 +
13.230 + excess.set(v, excess.get(v)+flo);
13.231 + exc-=flo;
13.232 + flow.set(e,0);
13.233 + }
13.234 + } else if ( newlevel > level.get(v) ) {
13.235 + newlevel = level.get(v);
13.236 + }
13.237 + } //for in edges vw
13.238 +
13.239 + } // if w still has excess after the out edge for cycle
13.240 +
13.241 + excess.set(w, exc);
13.242 +
13.243 + /*
13.244 + Relabel
13.245 + */
13.246 +
13.247 +
13.248 + if ( exc > 0 ) {
13.249 + //now 'lev' is the old level of w
13.250 +
13.251 + //unlacing starts
13.252 + NodeIt right_n=right.get(w);
13.253 + NodeIt left_n=left.get(w);
13.254 +
13.255 + if ( right_n != 0 ) {
13.256 + if ( left_n != 0 ) {
13.257 + right.set(left_n, right_n);
13.258 + left.set(right_n, left_n);
13.259 + } else {
13.260 + level_list[lev]=right_n;
13.261 + left.set(right_n, 0);
13.262 + }
13.263 + } else {
13.264 + if ( left_n != 0 ) {
13.265 + right.set(left_n, 0);
13.266 + } else {
13.267 + level_list[lev]=0;
13.268 +
13.269 + }
13.270 + }
13.271 + //unlacing ends
13.272 +
13.273 + //gapping starts
13.274 + if ( level_list[lev]==0 ) {
13.275 +
13.276 + for (int i=lev; i!=k ; ) {
13.277 + NodeIt v=level_list[++i];
13.278 + while ( v != 0 ) {
13.279 + level.set(v,n);
13.280 + v=right.get(v);
13.281 + }
13.282 + level_list[i]=0;
13.283 + if ( !what_heur ) active[i]=0;
13.284 + }
13.285 +
13.286 + level.set(w,n);
13.287 + b=lev-1;
13.288 + k=b;
13.289 + //gapping ends
13.290 + } else {
13.291 +
13.292 + if ( newlevel == n ) level.set(w,n);
13.293 + else {
13.294 + level.set(w,++newlevel);
13.295 + next.set(w,active[newlevel]);
13.296 + active[newlevel]=w;
13.297 + if ( what_heur ) b=newlevel;
13.298 + if ( k < newlevel ) ++k;
13.299 + NodeIt first=level_list[newlevel];
13.300 + if ( first != 0 ) left.set(first,w);
13.301 + right.set(w,first);
13.302 + left.set(w,0);
13.303 + level_list[newlevel]=w;
13.304 + }
13.305 + }
13.306 +
13.307 +
13.308 + ++relabel;
13.309 + if ( relabel >= heur ) {
13.310 + relabel=0;
13.311 + if ( what_heur ) {
13.312 + what_heur=0;
13.313 + heur=heur0;
13.314 + end=false;
13.315 + } else {
13.316 + what_heur=1;
13.317 + heur=heur1;
13.318 + b=k;
13.319 + }
13.320 + }
13.321 +
13.322 +
13.323 + } // if ( exc > 0 )
13.324 +
13.325 +
13.326 + } // if stack[b] is nonempty
13.327 +
13.328 + } // while(true)
13.329 +
13.330 +
13.331 +
13.332 + for( EachNodeIt v=G.template first<EachNodeIt>();
13.333 + v.valid(); ++v)
13.334 + if (level.get(v) >= n ) cut.set(v,true);
13.335 +
13.336 + value = excess.get(t);
13.337 + /*Max flow value.*/
13.338 +
13.339 + } //void run()
13.340 +
13.341 +
13.342 +
13.343 +
13.344 + T maxFlow() {
13.345 + return value;
13.346 + }
13.347 +
13.348 +
13.349 +
13.350 + CutMap minCut() {
13.351 + return cut;
13.352 + }
13.353 +
13.354 +
13.355 + };
13.356 +}//namespace
13.357 +#endif
13.358 +
13.359 +
13.360 +
13.361 +
14.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
14.2 +++ b/src/work/jacint/preflow_param.cc Sat Feb 21 21:01:22 2004 +0000
14.3 @@ -0,0 +1,47 @@
14.4 +#include <iostream>
14.5 +#include <fstream>
14.6 +
14.7 +#include <list_graph.hh>
14.8 +#include <dimacs.hh>
14.9 +#include <preflow_param.h>
14.10 +#include <time_measure.h>
14.11 +
14.12 +using namespace hugo;
14.13 +
14.14 +// Use a DIMACS max flow file as stdin.
14.15 +// read_dimacs_demo < dimacs_max_flow_file
14.16 +int main(int, char **) {
14.17 + typedef ListGraph::NodeIt NodeIt;
14.18 + typedef ListGraph::EachEdgeIt EachEdgeIt;
14.19 +
14.20 + ListGraph G;
14.21 + NodeIt s, t;
14.22 + ListGraph::EdgeMap<int> cap(G);
14.23 + readDimacsMaxFlow(std::cin, G, s, t, cap);
14.24 +
14.25 + std::cout << "preflow parameters test ..." << std::endl;
14.26 +
14.27 + double min=1000000;
14.28 +
14.29 + int _j=0;
14.30 + int _k=0;
14.31 +
14.32 + for (int j=1; j!=40; ++j ){
14.33 + for (int k=1; k!=j; ++k ){
14.34 +
14.35 + double mintime=1000000;
14.36 +
14.37 + for ( int i=1; i!=4; ++i ) {
14.38 + double pre_time=currTime();
14.39 + preflow_param<ListGraph, int> max_flow_test(G, s, t, cap, j, k);
14.40 + double post_time=currTime();
14.41 + if ( mintime > post_time-pre_time ) mintime = post_time-pre_time;
14.42 + }
14.43 + if (min > mintime ) {min=mintime; _j=j; _k=k;}
14.44 + }
14.45 + }
14.46 +
14.47 + std::cout<<"Min. at ("<<_j<<","<<_k<<") in time "<<min<<std::endl;
14.48 +
14.49 + return 0;
14.50 +}
15.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
15.2 +++ b/src/work/jacint/preflow_param.h Sat Feb 21 21:01:22 2004 +0000
15.3 @@ -0,0 +1,541 @@
15.4 +// -*- C++ -*-
15.5 +/*
15.6 +preflow_param.h
15.7 +by jacint.
15.8 +
15.9 +For testing the parameters H0, H1 of preflow.h
15.10 +
15.11 +The constructor runs the algorithm.
15.12 +
15.13 +Members:
15.14 +
15.15 +T maxFlow() : returns the value of a maximum flow
15.16 +
15.17 +T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
15.18 +
15.19 +FlowMap Flow() : returns the fixed maximum flow x
15.20 +
15.21 +void minMinCut(CutMap& M) : sets M to the characteristic vector of the
15.22 + minimum min cut. M should be a map of bools initialized to false.
15.23 +
15.24 +void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
15.25 + maximum min cut. M should be a map of bools initialized to false.
15.26 +
15.27 +void minCut(CutMap& M) : fast function, sets M to the characteristic
15.28 + vector of a minimum cut.
15.29 +
15.30 +Different member from the other preflow_hl-s (here we have a member
15.31 +'NodeMap<bool> cut').
15.32 +
15.33 +CutMap minCut() : fast function, giving the characteristic
15.34 + vector of a minimum cut.
15.35 +
15.36 +
15.37 +*/
15.38 +
15.39 +#ifndef PREFLOW_PARAM_H
15.40 +#define PREFLOW_PARAM_H
15.41 +
15.42 +#include <vector>
15.43 +#include <queue>
15.44 +
15.45 +#include <time_measure.h> //for test
15.46 +
15.47 +namespace hugo {
15.48 +
15.49 + template <typename Graph, typename T,
15.50 + typename FlowMap=typename Graph::EdgeMap<T>,
15.51 + typename CapMap=typename Graph::EdgeMap<T> >
15.52 + class preflow_param {
15.53 +
15.54 + typedef typename Graph::NodeIt NodeIt;
15.55 + typedef typename Graph::EdgeIt EdgeIt;
15.56 + typedef typename Graph::EachNodeIt EachNodeIt;
15.57 + typedef typename Graph::OutEdgeIt OutEdgeIt;
15.58 + typedef typename Graph::InEdgeIt InEdgeIt;
15.59 +
15.60 + Graph& G;
15.61 + NodeIt s;
15.62 + NodeIt t;
15.63 + FlowMap flow;
15.64 + CapMap& capacity;
15.65 + T value;
15.66 + int H0;
15.67 + int H1;
15.68 +
15.69 + public:
15.70 + double time;
15.71 +
15.72 + preflow_param(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity,
15.73 + int _H0, int _H1) :
15.74 + G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity), H0(_H0), H1(_H1) {
15.75 +
15.76 + bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd
15.77 + int n=G.nodeNum();
15.78 + int heur0=(int)(H0*n); //time while running 'bound decrease'
15.79 + int heur1=(int)(H1*n); //time while running 'highest label'
15.80 + int heur=heur1; //starting time interval (#of relabels)
15.81 + bool what_heur=1;
15.82 + /*
15.83 + what_heur is 0 in case 'bound decrease'
15.84 + and 1 in case 'highest label'
15.85 + */
15.86 + bool end=false; //for the 0 case
15.87 + /*
15.88 + Needed for 'bound decrease', 'true'
15.89 + means no active nodes are above bound b.
15.90 + */
15.91 + int relabel=0;
15.92 + int k=n-2;
15.93 + int b=k;
15.94 + /*
15.95 + b is a bound on the highest level of the stack.
15.96 + k is a bound on the highest nonempty level i < n.
15.97 + */
15.98 +
15.99 + typename Graph::NodeMap<int> level(G,n);
15.100 + typename Graph::NodeMap<T> excess(G);
15.101 +
15.102 + std::vector<NodeIt> active(n);
15.103 + typename Graph::NodeMap<NodeIt> next(G);
15.104 + //Stack of the active nodes in level i < n.
15.105 + //We use it in both phases.
15.106 +
15.107 + typename Graph::NodeMap<NodeIt> left(G);
15.108 + typename Graph::NodeMap<NodeIt> right(G);
15.109 + std::vector<NodeIt> level_list(n);
15.110 + /*
15.111 + Needed for the list of the nodes in level i.
15.112 + */
15.113 +
15.114 + /*Reverse_bfs from t, to find the starting level.*/
15.115 + level.set(t,0);
15.116 + std::queue<NodeIt> bfs_queue;
15.117 + bfs_queue.push(t);
15.118 +
15.119 + while (!bfs_queue.empty()) {
15.120 +
15.121 + NodeIt v=bfs_queue.front();
15.122 + bfs_queue.pop();
15.123 + int l=level.get(v)+1;
15.124 +
15.125 + for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
15.126 + NodeIt w=G.tail(e);
15.127 + if ( level.get(w) == n && w != s ) {
15.128 + bfs_queue.push(w);
15.129 + NodeIt first=level_list[l];
15.130 + if ( first != 0 ) left.set(first,w);
15.131 + right.set(w,first);
15.132 + level_list[l]=w;
15.133 + level.set(w, l);
15.134 + }
15.135 + }
15.136 + }
15.137 +
15.138 + level.set(s,n);
15.139 +
15.140 +
15.141 + /* Starting flow. It is everywhere 0 at the moment. */
15.142 + for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
15.143 + {
15.144 + T c=capacity.get(e);
15.145 + if ( c == 0 ) continue;
15.146 + NodeIt w=G.head(e);
15.147 + if ( level.get(w) < n ) {
15.148 + if ( excess.get(w) == 0 && w!=t ) {
15.149 + next.set(w,active[level.get(w)]);
15.150 + active[level.get(w)]=w;
15.151 + }
15.152 + flow.set(e, c);
15.153 + excess.set(w, excess.get(w)+c);
15.154 + }
15.155 + }
15.156 +
15.157 + /*
15.158 + End of preprocessing
15.159 + */
15.160 +
15.161 +
15.162 +
15.163 + /*
15.164 + Push/relabel on the highest level active nodes.
15.165 + */
15.166 + while ( true ) {
15.167 +
15.168 + if ( b == 0 ) {
15.169 + if ( phase ) break;
15.170 +
15.171 + if ( !what_heur && !end && k > 0 ) {
15.172 + b=k;
15.173 + end=true;
15.174 + } else {
15.175 + time=currTime();
15.176 + phase=1;
15.177 + level.set(s,0);
15.178 +
15.179 + std::queue<NodeIt> bfs_queue;
15.180 + bfs_queue.push(s);
15.181 +
15.182 + while (!bfs_queue.empty()) {
15.183 +
15.184 + NodeIt v=bfs_queue.front();
15.185 + bfs_queue.pop();
15.186 + int l=level.get(v)+1;
15.187 +
15.188 + for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
15.189 + if ( capacity.get(e) == flow.get(e) ) continue;
15.190 + NodeIt u=G.tail(e);
15.191 + if ( level.get(u) >= n ) {
15.192 + bfs_queue.push(u);
15.193 + level.set(u, l);
15.194 + if ( excess.get(u) > 0 ) {
15.195 + next.set(u,active[l]);
15.196 + active[l]=u;
15.197 + }
15.198 + }
15.199 + }
15.200 +
15.201 +
15.202 + for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
15.203 + if ( 0 == flow.get(e) ) continue;
15.204 + NodeIt u=G.head(e);
15.205 + if ( level.get(u) >= n ) {
15.206 + bfs_queue.push(u);
15.207 + level.set(u, l);
15.208 + if ( excess.get(u) > 0 ) {
15.209 + next.set(u,active[l]);
15.210 + active[l]=u;
15.211 + }
15.212 + }
15.213 + }
15.214 + }
15.215 + b=n-2;
15.216 + }
15.217 +
15.218 + }
15.219 +
15.220 +
15.221 + if ( active[b] == 0 ) --b;
15.222 + else {
15.223 + end=false; //needed only for phase 0, case hl2
15.224 +
15.225 + NodeIt w=active[b]; //w is a highest label active node.
15.226 + active[b]=next.get(w);
15.227 + int lev=level.get(w);
15.228 + T exc=excess.get(w);
15.229 + int newlevel=n; //In newlevel we bound the next level of w.
15.230 +
15.231 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
15.232 +
15.233 + if ( flow.get(e) == capacity.get(e) ) continue;
15.234 + NodeIt v=G.head(e);
15.235 + //e=wv
15.236 +
15.237 + if( lev > level.get(v) ) {
15.238 + /*Push is allowed now*/
15.239 +
15.240 + if ( excess.get(v)==0 && v!=t && v!=s ) {
15.241 + int lev_v=level.get(v);
15.242 + next.set(v,active[lev_v]);
15.243 + active[lev_v]=v;
15.244 + }
15.245 + /*v becomes active.*/
15.246 +
15.247 + T cap=capacity.get(e);
15.248 + T flo=flow.get(e);
15.249 + T remcap=cap-flo;
15.250 +
15.251 + if ( remcap >= exc ) {
15.252 + /*A nonsaturating push.*/
15.253 +
15.254 + flow.set(e, flo+exc);
15.255 + excess.set(v, excess.get(v)+exc);
15.256 + exc=0;
15.257 + break;
15.258 +
15.259 + } else {
15.260 + /*A saturating push.*/
15.261 +
15.262 + flow.set(e, cap);
15.263 + excess.set(v, excess.get(v)+remcap);
15.264 + exc-=remcap;
15.265 + }
15.266 + } else if ( newlevel > level.get(v) ){
15.267 + newlevel = level.get(v);
15.268 + }
15.269 +
15.270 + } //for out edges wv
15.271 +
15.272 +
15.273 + if ( exc > 0 ) {
15.274 + for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
15.275 +
15.276 + if( flow.get(e) == 0 ) continue;
15.277 + NodeIt v=G.tail(e);
15.278 + //e=vw
15.279 +
15.280 + if( lev > level.get(v) ) {
15.281 + /*Push is allowed now*/
15.282 +
15.283 + if ( excess.get(v)==0 && v!=t && v!=s ) {
15.284 + int lev_v=level.get(v);
15.285 + next.set(v,active[lev_v]);
15.286 + active[lev_v]=v;
15.287 + /*v becomes active.*/
15.288 + }
15.289 +
15.290 + T flo=flow.get(e);
15.291 +
15.292 + if ( flo >= exc ) {
15.293 + /*A nonsaturating push.*/
15.294 +
15.295 + flow.set(e, flo-exc);
15.296 + excess.set(v, excess.get(v)+exc);
15.297 + exc=0;
15.298 + break;
15.299 + } else {
15.300 + /*A saturating push.*/
15.301 +
15.302 + excess.set(v, excess.get(v)+flo);
15.303 + exc-=flo;
15.304 + flow.set(e,0);
15.305 + }
15.306 + } else if ( newlevel > level.get(v) ) {
15.307 + newlevel = level.get(v);
15.308 + }
15.309 + } //for in edges vw
15.310 +
15.311 + } // if w still has excess after the out edge for cycle
15.312 +
15.313 + excess.set(w, exc);
15.314 +
15.315 + /*
15.316 + Relabel
15.317 + */
15.318 +
15.319 +
15.320 + if ( exc > 0 ) {
15.321 + //now 'lev' is the old level of w
15.322 +
15.323 + if ( phase ) {
15.324 + level.set(w,++newlevel);
15.325 + next.set(w,active[newlevel]);
15.326 + active[newlevel]=w;
15.327 + b=newlevel;
15.328 + } else {
15.329 + //unlacing
15.330 + NodeIt right_n=right.get(w);
15.331 + NodeIt left_n=left.get(w);
15.332 +
15.333 + if ( right_n != 0 ) {
15.334 + if ( left_n != 0 ) {
15.335 + right.set(left_n, right_n);
15.336 + left.set(right_n, left_n);
15.337 + } else {
15.338 + level_list[lev]=right_n;
15.339 + left.set(right_n, 0);
15.340 + }
15.341 + } else {
15.342 + if ( left_n != 0 ) {
15.343 + right.set(left_n, 0);
15.344 + } else {
15.345 + level_list[lev]=0;
15.346 +
15.347 + }
15.348 + }
15.349 +
15.350 + if ( level_list[lev]==0 ) {
15.351 +
15.352 + for (int i=lev; i!=k ; ) {
15.353 + NodeIt v=level_list[++i];
15.354 + while ( v != 0 ) {
15.355 + level.set(v,n);
15.356 + v=right.get(v);
15.357 + }
15.358 + level_list[i]=0;
15.359 + if ( !what_heur ) active[i]=0;
15.360 + }
15.361 +
15.362 + level.set(w,n);
15.363 + b=lev-1;
15.364 + k=b;
15.365 +
15.366 + } else {
15.367 +
15.368 + if ( newlevel == n ) level.set(w,n);
15.369 + else {
15.370 + level.set(w,++newlevel);
15.371 + next.set(w,active[newlevel]);
15.372 + active[newlevel]=w;
15.373 + if ( what_heur ) b=newlevel;
15.374 + if ( k < newlevel ) ++k;
15.375 + NodeIt first=level_list[newlevel];
15.376 + if ( first != 0 ) left.set(first,w);
15.377 + right.set(w,first);
15.378 + left.set(w,0);
15.379 + level_list[newlevel]=w;
15.380 + }
15.381 + }
15.382 +
15.383 +
15.384 + ++relabel;
15.385 + if ( relabel >= heur ) {
15.386 + relabel=0;
15.387 + if ( what_heur ) {
15.388 + what_heur=0;
15.389 + heur=heur0;
15.390 + end=false;
15.391 + } else {
15.392 + what_heur=1;
15.393 + heur=heur1;
15.394 + b=k;
15.395 + }
15.396 + }
15.397 + } //phase 0
15.398 +
15.399 +
15.400 + } // if ( exc > 0 )
15.401 +
15.402 +
15.403 + } // if stack[b] is nonempty
15.404 +
15.405 + } // while(true)
15.406 +
15.407 +
15.408 + value = excess.get(t);
15.409 + /*Max flow value.*/
15.410 +
15.411 + } //void run()
15.412 +
15.413 +
15.414 +
15.415 +
15.416 +
15.417 + /*
15.418 + Returns the maximum value of a flow.
15.419 + */
15.420 +
15.421 + T maxFlow() {
15.422 + return value;
15.423 + }
15.424 +
15.425 +
15.426 +
15.427 + /*
15.428 + For the maximum flow x found by the algorithm,
15.429 + it returns the flow value on edge e, i.e. x(e).
15.430 + */
15.431 +
15.432 + T flowOnEdge(EdgeIt e) {
15.433 + return flow.get(e);
15.434 + }
15.435 +
15.436 +
15.437 +
15.438 + FlowMap Flow() {
15.439 + return flow;
15.440 + }
15.441 +
15.442 +
15.443 +
15.444 + void Flow(FlowMap& _flow ) {
15.445 + for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v)
15.446 + _flow.set(v,flow.get(v));
15.447 + }
15.448 +
15.449 +
15.450 +
15.451 + /*
15.452 + Returns the minimum min cut, by a bfs from s in the residual graph.
15.453 + */
15.454 +
15.455 + template<typename CutMap>
15.456 + void minCut(CutMap& M) {
15.457 +
15.458 + std::queue<NodeIt> queue;
15.459 +
15.460 + M.set(s,true);
15.461 + queue.push(s);
15.462 +
15.463 + while (!queue.empty()) {
15.464 + NodeIt w=queue.front();
15.465 + queue.pop();
15.466 +
15.467 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
15.468 + NodeIt v=G.head(e);
15.469 + if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
15.470 + queue.push(v);
15.471 + M.set(v, true);
15.472 + }
15.473 + }
15.474 +
15.475 + for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
15.476 + NodeIt v=G.tail(e);
15.477 + if (!M.get(v) && flow.get(e) > 0 ) {
15.478 + queue.push(v);
15.479 + M.set(v, true);
15.480 + }
15.481 + }
15.482 +
15.483 + }
15.484 +
15.485 + }
15.486 +
15.487 +
15.488 +
15.489 + /*
15.490 + Returns the maximum min cut, by a reverse bfs
15.491 + from t in the residual graph.
15.492 + */
15.493 +
15.494 + template<typename CutMap>
15.495 + void maxMinCut(CutMap& M) {
15.496 +
15.497 + std::queue<NodeIt> queue;
15.498 +
15.499 + M.set(t,true);
15.500 + queue.push(t);
15.501 +
15.502 + while (!queue.empty()) {
15.503 + NodeIt w=queue.front();
15.504 + queue.pop();
15.505 +
15.506 + for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
15.507 + NodeIt v=G.tail(e);
15.508 + if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
15.509 + queue.push(v);
15.510 + M.set(v, true);
15.511 + }
15.512 + }
15.513 +
15.514 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
15.515 + NodeIt v=G.head(e);
15.516 + if (!M.get(v) && flow.get(e) > 0 ) {
15.517 + queue.push(v);
15.518 + M.set(v, true);
15.519 + }
15.520 + }
15.521 + }
15.522 +
15.523 + for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) {
15.524 + M.set(v, !M.get(v));
15.525 + }
15.526 +
15.527 + }
15.528 +
15.529 +
15.530 +
15.531 + template<typename CutMap>
15.532 + void minMinCut(CutMap& M) {
15.533 + minCut(M);
15.534 + }
15.535 +
15.536 +
15.537 +
15.538 + };
15.539 +}//namespace
15.540 +#endif
15.541 +
15.542 +
15.543 +
15.544 +
16.1 --- a/src/work/jacint/preflow_push_hl.h Fri Feb 20 22:01:02 2004 +0000
16.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
16.3 @@ -1,398 +0,0 @@
16.4 -// -*- C++ -*-
16.5 -
16.6 -//kerdesek: nem tudom lehet-e a
16.7 -//kieleket csak a legf n szintu pontokra nezni.
16.8 -
16.9 -/*
16.10 -preflow_push_hl.h
16.11 -by jacint.
16.12 -Runs the highest label variant of the preflow push algorithm with
16.13 -running time O(n^2\sqrt(m)), and with the 'empty level' heuristic.
16.14 -
16.15 -'A' is a parameter for the empty_level heuristic
16.16 -
16.17 -Member functions:
16.18 -
16.19 -void run() : runs the algorithm
16.20 -
16.21 - The following functions should be used after run() was already run.
16.22 -
16.23 -T maxflow() : returns the value of a maximum flow
16.24 -
16.25 -T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
16.26 -
16.27 -FlowMap allflow() : returns the fixed maximum flow x
16.28 -
16.29 -void mincut(CutMap& M) : sets M to the characteristic vector of a
16.30 - minimum cut. M should be a map of bools initialized to false.
16.31 -
16.32 -void min_mincut(CutMap& M) : sets M to the characteristic vector of the
16.33 - minimum min cut. M should be a map of bools initialized to false.
16.34 -
16.35 -void max_mincut(CutMap& M) : sets M to the characteristic vector of the
16.36 - maximum min cut. M should be a map of bools initialized to false.
16.37 -
16.38 -*/
16.39 -
16.40 -#ifndef PREFLOW_PUSH_HL_H
16.41 -#define PREFLOW_PUSH_HL_H
16.42 -
16.43 -#define A 1
16.44 -
16.45 -#include <vector>
16.46 -#include <stack>
16.47 -#include <queue>
16.48 -
16.49 -namespace hugo {
16.50 -
16.51 - template <typename Graph, typename T,
16.52 - typename FlowMap=typename Graph::EdgeMap<T>, typename CapMap=typename Graph::EdgeMap<T>,
16.53 - typename IntMap=typename Graph::NodeMap<int>, typename TMap=typename Graph::NodeMap<T> >
16.54 - class preflow_push_hl {
16.55 -
16.56 - typedef typename Graph::NodeIt NodeIt;
16.57 - typedef typename Graph::EdgeIt EdgeIt;
16.58 - typedef typename Graph::EachNodeIt EachNodeIt;
16.59 - typedef typename Graph::OutEdgeIt OutEdgeIt;
16.60 - typedef typename Graph::InEdgeIt InEdgeIt;
16.61 -
16.62 - Graph& G;
16.63 - NodeIt s;
16.64 - NodeIt t;
16.65 - FlowMap flow;
16.66 - CapMap& capacity;
16.67 - T value;
16.68 -
16.69 - public:
16.70 -
16.71 - preflow_push_hl(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
16.72 - G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) { }
16.73 -
16.74 -
16.75 -
16.76 -
16.77 - void run() {
16.78 -
16.79 - int n=G.nodeNum();
16.80 - int b=n-2;
16.81 - /*
16.82 - b is a bound on the highest level of an active node.
16.83 - */
16.84 -
16.85 - IntMap level(G,n);
16.86 - TMap excess(G);
16.87 -
16.88 - std::vector<int> numb(n);
16.89 - /*
16.90 - The number of nodes on level i < n. It is
16.91 - initialized to n+1, because of the reverse_bfs-part.
16.92 - */
16.93 -
16.94 - std::vector<std::stack<NodeIt> > stack(2*n-1);
16.95 - //Stack of the active nodes in level i.
16.96 -
16.97 -
16.98 - /*Reverse_bfs from t, to find the starting level.*/
16.99 - level.set(t,0);
16.100 - std::queue<NodeIt> bfs_queue;
16.101 - bfs_queue.push(t);
16.102 -
16.103 - while (!bfs_queue.empty()) {
16.104 -
16.105 - NodeIt v=bfs_queue.front();
16.106 - bfs_queue.pop();
16.107 - int l=level.get(v)+1;
16.108 -
16.109 - for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
16.110 - NodeIt w=G.tail(e);
16.111 - if ( level.get(w) == n ) {
16.112 - bfs_queue.push(w);
16.113 - ++numb[l];
16.114 - level.set(w, l);
16.115 - }
16.116 - }
16.117 - }
16.118 -
16.119 - level.set(s,n);
16.120 -
16.121 -
16.122 -
16.123 - /* Starting flow. It is everywhere 0 at the moment. */
16.124 - for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
16.125 - {
16.126 - T c=capacity.get(e);
16.127 - if ( c == 0 ) continue;
16.128 - NodeIt w=G.head(e);
16.129 - if ( w!=s ) {
16.130 - if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w);
16.131 - flow.set(e, c);
16.132 - excess.set(w, excess.get(w)+c);
16.133 - }
16.134 - }
16.135 -
16.136 - /*
16.137 - End of preprocessing
16.138 - */
16.139 -
16.140 -
16.141 - /*
16.142 - Push/relabel on the highest level active nodes.
16.143 - */
16.144 - /*While there exists an active node.*/
16.145 - while (b) {
16.146 - if ( stack[b].empty() ) {
16.147 - --b;
16.148 - continue;
16.149 - }
16.150 -
16.151 - NodeIt w=stack[b].top(); //w is a highest label active node.
16.152 - stack[b].pop();
16.153 - int lev=level.get(w);
16.154 - int exc=excess.get(w);
16.155 - int newlevel=2*n; //In newlevel we bound the next level of w.
16.156 - //vagy MAXINT
16.157 -
16.158 - // if ( level.get(w) < n ) { //Nem tudom ez mukodik-e
16.159 - for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
16.160 -
16.161 - if ( flow.get(e) == capacity.get(e) ) continue;
16.162 - NodeIt v=G.head(e);
16.163 - //e=wv
16.164 -
16.165 - if( lev > level.get(v) ) {
16.166 - /*Push is allowed now*/
16.167 -
16.168 - if ( excess.get(v)==0 && v != s && v !=t )
16.169 - stack[level.get(v)].push(v);
16.170 - /*v becomes active.*/
16.171 -
16.172 - int cap=capacity.get(e);
16.173 - int flo=flow.get(e);
16.174 - int remcap=cap-flo;
16.175 -
16.176 - if ( remcap >= exc ) {
16.177 - /*A nonsaturating push.*/
16.178 -
16.179 - flow.set(e, flo+exc);
16.180 - excess.set(v, excess.get(v)+exc);
16.181 - exc=0;
16.182 - break;
16.183 -
16.184 - } else {
16.185 - /*A saturating push.*/
16.186 -
16.187 - flow.set(e, cap );
16.188 - excess.set(v, excess.get(v)+remcap);
16.189 - exc-=remcap;
16.190 - }
16.191 - } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
16.192 -
16.193 - } //for out edges wv
16.194 -
16.195 -
16.196 - if ( exc > 0 ) {
16.197 - for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
16.198 -
16.199 - if( flow.get(e) == 0 ) continue;
16.200 - NodeIt v=G.tail(e);
16.201 - //e=vw
16.202 -
16.203 - if( lev > level.get(v) ) {
16.204 - /*Push is allowed now*/
16.205 -
16.206 - if ( excess.get(v)==0 && v != s && v !=t)
16.207 - stack[level.get(v)].push(v);
16.208 - /*v becomes active.*/
16.209 -
16.210 - int flo=flow.get(e);
16.211 -
16.212 - if ( flo >= exc ) {
16.213 - /*A nonsaturating push.*/
16.214 -
16.215 - flow.set(e, flo-exc);
16.216 - excess.set(v, excess.get(v)+exc);
16.217 - exc=0;
16.218 - break;
16.219 - } else {
16.220 - /*A saturating push.*/
16.221 -
16.222 - excess.set(v, excess.get(v)+flo);
16.223 - exc-=flo;
16.224 - flow.set(e,0);
16.225 - }
16.226 - } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
16.227 -
16.228 - } //for in edges vw
16.229 -
16.230 - } // if w still has excess after the out edge for cycle
16.231 -
16.232 - excess.set(w, exc);
16.233 -
16.234 -
16.235 -
16.236 -
16.237 - /*
16.238 - Relabel
16.239 - */
16.240 -
16.241 - if ( exc > 0 ) {
16.242 - //now 'lev' is the old level of w
16.243 - level.set(w,++newlevel);
16.244 -
16.245 - if ( lev < n ) {
16.246 - --numb[lev];
16.247 -
16.248 - if ( !numb[lev] && lev < A*n ) { //If the level of w gets empty.
16.249 -
16.250 - for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
16.251 - if (level.get(v) > lev && level.get(v) < n ) level.set(v,n);
16.252 - }
16.253 - for (int i=lev+1 ; i!=n ; ++i) numb[i]=0;
16.254 - if ( newlevel < n ) newlevel=n;
16.255 - } else {
16.256 - if ( newlevel < n ) ++numb[newlevel];
16.257 - }
16.258 - }
16.259 -
16.260 - stack[newlevel].push(w);
16.261 - b=newlevel;
16.262 -
16.263 - }
16.264 -
16.265 - } // while(b)
16.266 -
16.267 -
16.268 - value = excess.get(t);
16.269 - /*Max flow value.*/
16.270 -
16.271 -
16.272 - } //void run()
16.273 -
16.274 -
16.275 -
16.276 -
16.277 -
16.278 - /*
16.279 - Returns the maximum value of a flow.
16.280 - */
16.281 -
16.282 - T maxflow() {
16.283 - return value;
16.284 - }
16.285 -
16.286 -
16.287 -
16.288 - /*
16.289 - For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e).
16.290 - */
16.291 -
16.292 - T flowonedge(const EdgeIt e) {
16.293 - return flow.get(e);
16.294 - }
16.295 -
16.296 -
16.297 -
16.298 - /*
16.299 - Returns the maximum flow x found by the algorithm.
16.300 - */
16.301 -
16.302 - FlowMap allflow() {
16.303 - return flow;
16.304 - }
16.305 -
16.306 -
16.307 -
16.308 -
16.309 - /*
16.310 - Returns the minimum min cut, by a bfs from s in the residual graph.
16.311 - */
16.312 -
16.313 - template<typename CutMap>
16.314 - void mincut(CutMap& M) {
16.315 -
16.316 - std::queue<NodeIt> queue;
16.317 -
16.318 - M.set(s,true);
16.319 - queue.push(s);
16.320 -
16.321 - while (!queue.empty()) {
16.322 - NodeIt w=queue.front();
16.323 - queue.pop();
16.324 -
16.325 - for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
16.326 - NodeIt v=G.head(e);
16.327 - if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
16.328 - queue.push(v);
16.329 - M.set(v, true);
16.330 - }
16.331 - }
16.332 -
16.333 - for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
16.334 - NodeIt v=G.tail(e);
16.335 - if (!M.get(v) && flow.get(e) > 0 ) {
16.336 - queue.push(v);
16.337 - M.set(v, true);
16.338 - }
16.339 - }
16.340 -
16.341 - }
16.342 - }
16.343 -
16.344 -
16.345 -
16.346 - /*
16.347 - Returns the maximum min cut, by a reverse bfs
16.348 - from t in the residual graph.
16.349 - */
16.350 -
16.351 - template<typename CutMap>
16.352 - void max_mincut(CutMap& M) {
16.353 -
16.354 - std::queue<NodeIt> queue;
16.355 -
16.356 - M.set(t,true);
16.357 - queue.push(t);
16.358 -
16.359 - while (!queue.empty()) {
16.360 - NodeIt w=queue.front();
16.361 - queue.pop();
16.362 -
16.363 - for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
16.364 - NodeIt v=G.tail(e);
16.365 - if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
16.366 - queue.push(v);
16.367 - M.set(v, true);
16.368 - }
16.369 - }
16.370 -
16.371 - for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
16.372 - NodeIt v=G.head(e);
16.373 - if (!M.get(v) && flow.get(e) > 0 ) {
16.374 - queue.push(v);
16.375 - M.set(v, true);
16.376 - }
16.377 - }
16.378 - }
16.379 -
16.380 - for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) {
16.381 - M.set(v, !M.get(v));
16.382 - }
16.383 -
16.384 - }
16.385 -
16.386 -
16.387 -
16.388 - template<typename CutMap>
16.389 - void min_mincut(CutMap& M) {
16.390 - mincut(M);
16.391 - }
16.392 -
16.393 -
16.394 -
16.395 - };
16.396 -}//namespace hugo
16.397 -#endif
16.398 -
16.399 -
16.400 -
16.401 -
17.1 --- a/src/work/jacint/preflow_push_hl.hh Fri Feb 20 22:01:02 2004 +0000
17.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
17.3 @@ -1,320 +0,0 @@
17.4 -/*
17.5 -preflow_push_hl.hh
17.6 -by jacint.
17.7 -Runs the highest label variant of the preflow push algorithm with
17.8 -running time O(n^2\sqrt(m)).
17.9 -
17.10 -Member functions:
17.11 -
17.12 -void run() : runs the algorithm
17.13 -
17.14 - The following functions should be used after run() was already run.
17.15 -
17.16 -T maxflow() : returns the value of a maximum flow
17.17 -
17.18 -T flowonedge(edge_iterator e) : for a fixed maximum flow x it returns x(e)
17.19 -
17.20 -edge_property_vector<graph_type, T> allflow() : returns the fixed maximum flow x
17.21 -
17.22 -node_property_vector<graph_type, bool> mincut() : returns a
17.23 - characteristic vector of a minimum cut. (An empty level
17.24 - in the algorithm gives a minimum cut.)
17.25 -*/
17.26 -
17.27 -#ifndef PREFLOW_PUSH_HL_HH
17.28 -#define PREFLOW_PUSH_HL_HH
17.29 -
17.30 -#include <algorithm>
17.31 -#include <vector>
17.32 -#include <stack>
17.33 -
17.34 -#include <marci_graph_traits.hh>
17.35 -#include <marci_property_vector.hh>
17.36 -#include <reverse_bfs.hh>
17.37 -
17.38 -namespace hugo {
17.39 -
17.40 - template <typename graph_type, typename T>
17.41 - class preflow_push_hl {
17.42 -
17.43 - typedef typename graph_traits<graph_type>::node_iterator node_iterator;
17.44 - typedef typename graph_traits<graph_type>::edge_iterator edge_iterator;
17.45 - typedef typename graph_traits<graph_type>::each_node_iterator each_node_iterator;
17.46 - typedef typename graph_traits<graph_type>::out_edge_iterator out_edge_iterator;
17.47 - typedef typename graph_traits<graph_type>::in_edge_iterator in_edge_iterator;
17.48 - typedef typename graph_traits<graph_type>::each_edge_iterator each_edge_iterator;
17.49 -
17.50 -
17.51 - graph_type& G;
17.52 - node_iterator s;
17.53 - node_iterator t;
17.54 - edge_property_vector<graph_type, T> flow;
17.55 - edge_property_vector<graph_type, T>& capacity;
17.56 - T value;
17.57 - node_property_vector<graph_type, bool> mincutvector;
17.58 -
17.59 -
17.60 - public:
17.61 -
17.62 - preflow_push_hl(graph_type& _G, node_iterator _s, node_iterator _t, edge_property_vector<graph_type, T>& _capacity) : G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity), mincutvector(_G, true) { }
17.63 -
17.64 -
17.65 -
17.66 -
17.67 - /*
17.68 - The run() function runs the highest label preflow-push,
17.69 - running time: O(n^2\sqrt(m))
17.70 - */
17.71 - void run() {
17.72 -
17.73 - node_property_vector<graph_type, int> level(G); //level of node
17.74 - node_property_vector<graph_type, T> excess(G); //excess of node
17.75 -
17.76 - int n=number_of(G.first_node()); //number of nodes
17.77 - int b=n;
17.78 - /*b is a bound on the highest level of an active node. In the beginning it is at most n-2.*/
17.79 -
17.80 - std::vector<std::stack<node_iterator> > stack(2*n-1); //Stack of the active nodes in level i.
17.81 -
17.82 -
17.83 -
17.84 -
17.85 - /*Reverse_bfs from t, to find the starting level.*/
17.86 -
17.87 - reverse_bfs<list_graph> bfs(G, t);
17.88 - bfs.run();
17.89 - for(each_node_iterator v=G.first_node(); v.valid(); ++v) {
17.90 - level.put(v, bfs.dist(v));
17.91 - //std::cout << "the level of " << v << " is " << bfs.dist(v);
17.92 - }
17.93 -
17.94 - /*The level of s is fixed to n*/
17.95 - level.put(s,n);
17.96 -
17.97 -
17.98 -
17.99 -
17.100 -
17.101 - /* Starting flow. It is everywhere 0 at the moment. */
17.102 -
17.103 - for(out_edge_iterator i=G.first_out_edge(s); i.valid(); ++i)
17.104 - {
17.105 - node_iterator w=G.head(i);
17.106 - flow.put(i, capacity.get(i));
17.107 - stack[bfs.dist(w)].push(w);
17.108 - excess.put(w, capacity.get(i));
17.109 - }
17.110 -
17.111 -
17.112 - /*
17.113 - End of preprocessing
17.114 - */
17.115 -
17.116 -
17.117 -
17.118 - /*
17.119 - Push/relabel on the highest level active nodes.
17.120 - */
17.121 -
17.122 - /*While there exists active node.*/
17.123 - while (b) {
17.124 -
17.125 - /*We decrease the bound if there is no active node of level b.*/
17.126 - if (stack[b].empty()) {
17.127 - --b;
17.128 - } else {
17.129 -
17.130 - node_iterator w=stack[b].top(); //w is the highest label active node.
17.131 - stack[b].pop(); //We delete w from the stack.
17.132 -
17.133 - int newlevel=2*n-2; //In newlevel we maintain the next level of w.
17.134 -
17.135 - for(out_edge_iterator e=G.first_out_edge(w); e.valid(); ++e) {
17.136 - node_iterator v=G.head(e);
17.137 - /*e is the edge wv.*/
17.138 -
17.139 - if (flow.get(e)<capacity.get(e)) {
17.140 - /*e is an edge of the residual graph */
17.141 -
17.142 - if(level.get(w)==level.get(v)+1) {
17.143 - /*Push is allowed now*/
17.144 -
17.145 - if (capacity.get(e)-flow.get(e) > excess.get(w)) {
17.146 - /*A nonsaturating push.*/
17.147 -
17.148 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
17.149 - /*v becomes active.*/
17.150 -
17.151 - flow.put(e, flow.get(e)+excess.get(w));
17.152 - excess.put(v, excess.get(v)+excess.get(w));
17.153 - excess.put(w,0);
17.154 - //std::cout << w << " " << v <<" elore elen nonsat pump " << std::endl;
17.155 - break;
17.156 - } else {
17.157 - /*A saturating push.*/
17.158 -
17.159 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
17.160 - /*v becomes active.*/
17.161 -
17.162 - excess.put(v, excess.get(v)+capacity.get(e)-flow.get(e));
17.163 - excess.put(w, excess.get(w)-capacity.get(e)+flow.get(e));
17.164 - flow.put(e, capacity.get(e));
17.165 - //std::cout << w<<" " <<v<<" elore elen sat pump " << std::endl;
17.166 - if (excess.get(w)==0) break;
17.167 - /*If w is not active any more, then we go on to the next node.*/
17.168 -
17.169 - } // if (capacity.get(e)-flow.get(e) > excess.get(w))
17.170 - } // if(level.get(w)==level.get(v)+1)
17.171 -
17.172 - else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);}
17.173 -
17.174 - } //if (flow.get(e)<capacity.get(e))
17.175 -
17.176 - } //for(out_edge_iterator e=G.first_out_edge(w); e.valid(); ++e)
17.177 -
17.178 -
17.179 -
17.180 - for(in_edge_iterator e=G.first_in_edge(w); e.valid(); ++e) {
17.181 - node_iterator v=G.tail(e);
17.182 - /*e is the edge vw.*/
17.183 -
17.184 - if (excess.get(w)==0) break;
17.185 - /*It may happen, that w became inactive in the first for cycle.*/
17.186 - if(flow.get(e)>0) {
17.187 - /*e is an edge of the residual graph */
17.188 -
17.189 - if(level.get(w)==level.get(v)+1) {
17.190 - /*Push is allowed now*/
17.191 -
17.192 - if (flow.get(e) > excess.get(w)) {
17.193 - /*A nonsaturating push.*/
17.194 -
17.195 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
17.196 - /*v becomes active.*/
17.197 -
17.198 - flow.put(e, flow.get(e)-excess.get(w));
17.199 - excess.put(v, excess.get(v)+excess.get(w));
17.200 - excess.put(w,0);
17.201 - //std::cout << v << " " << w << " vissza elen nonsat pump " << std::endl;
17.202 - break;
17.203 - } else {
17.204 - /*A saturating push.*/
17.205 -
17.206 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
17.207 - /*v becomes active.*/
17.208 -
17.209 - excess.put(v, excess.get(v)+flow.get(e));
17.210 - excess.put(w, excess.get(w)-flow.get(e));
17.211 - flow.put(e,0);
17.212 - //std::cout << v <<" " << w << " vissza elen sat pump " << std::endl;
17.213 - if (excess.get(w)==0) { break;}
17.214 - } //if (flow.get(e) > excess.get(v))
17.215 - } //if(level.get(w)==level.get(v)+1)
17.216 -
17.217 - else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);}
17.218 -
17.219 -
17.220 - } //if (flow.get(e)>0)
17.221 -
17.222 - } //for
17.223 -
17.224 -
17.225 - if (excess.get(w)>0) {
17.226 - level.put(w,++newlevel);
17.227 - stack[newlevel].push(w);
17.228 - b=newlevel;
17.229 - //std::cout << "The new level of " << w << " is "<< newlevel <<std::endl;
17.230 - }
17.231 -
17.232 -
17.233 - } //else
17.234 -
17.235 - } //while(b)
17.236 -
17.237 - value = excess.get(t);
17.238 - /*Max flow value.*/
17.239 -
17.240 -
17.241 -
17.242 -
17.243 - } //void run()
17.244 -
17.245 -
17.246 -
17.247 -
17.248 -
17.249 - /*
17.250 - Returns the maximum value of a flow.
17.251 - */
17.252 -
17.253 - T maxflow() {
17.254 - return value;
17.255 - }
17.256 -
17.257 -
17.258 -
17.259 - /*
17.260 - For the maximum flow x found by the algorithm, it returns the flow value on edge e, i.e. x(e).
17.261 - */
17.262 -
17.263 - T flowonedge(edge_iterator e) {
17.264 - return flow.get(e);
17.265 - }
17.266 -
17.267 -
17.268 -
17.269 - /*
17.270 - Returns the maximum flow x found by the algorithm.
17.271 - */
17.272 -
17.273 - edge_property_vector<graph_type, T> allflow() {
17.274 - return flow;
17.275 - }
17.276 -
17.277 -
17.278 -
17.279 - /*
17.280 - Returns a minimum cut by using a reverse bfs from t in the residual graph.
17.281 - */
17.282 -
17.283 - node_property_vector<graph_type, bool> mincut() {
17.284 -
17.285 - std::queue<node_iterator> queue;
17.286 -
17.287 - mincutvector.put(t,false);
17.288 - queue.push(t);
17.289 -
17.290 - while (!queue.empty()) {
17.291 - node_iterator w=queue.front();
17.292 - queue.pop();
17.293 -
17.294 - for(in_edge_iterator e=G.first_in_edge(w) ; e.valid(); ++e) {
17.295 - node_iterator v=G.tail(e);
17.296 - if (mincutvector.get(v) && flow.get(e) < capacity.get(e) ) {
17.297 - queue.push(v);
17.298 - mincutvector.put(v, false);
17.299 - }
17.300 - } // for
17.301 -
17.302 - for(out_edge_iterator e=G.first_out_edge(w) ; e.valid(); ++e) {
17.303 - node_iterator v=G.head(e);
17.304 - if (mincutvector.get(v) && flow.get(e) > 0 ) {
17.305 - queue.push(v);
17.306 - mincutvector.put(v, false);
17.307 - }
17.308 - } // for
17.309 -
17.310 - }
17.311 -
17.312 - return mincutvector;
17.313 -
17.314 - }
17.315 -
17.316 -
17.317 - };
17.318 -}//namespace hugo
17.319 -#endif
17.320 -
17.321 -
17.322 -
17.323 -
18.1 --- a/src/work/jacint/preflow_push_max_flow.h Fri Feb 20 22:01:02 2004 +0000
18.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
18.3 @@ -1,296 +0,0 @@
18.4 -// -*- C++ -*-
18.5 -/*
18.6 -preflow_push_max_flow_h
18.7 -by jacint.
18.8 -Runs a preflow push algorithm with the modification,
18.9 -that we do not push on nodes with level at least n.
18.10 -Moreover, if a level gets empty, we set all nodes above that
18.11 -level to level n. Hence, in the end, we arrive at a maximum preflow
18.12 -with value of a max flow value. An empty level gives a minimum cut.
18.13 -
18.14 -Member functions:
18.15 -
18.16 -void run() : runs the algorithm
18.17 -
18.18 - The following functions should be used after run() was already run.
18.19 -
18.20 -T maxflow() : returns the value of a maximum flow
18.21 -
18.22 -void mincut(CutMap& M) : sets M to the characteristic vector of a
18.23 - minimum cut. M should be a map of bools initialized to false.
18.24 -
18.25 -*/
18.26 -
18.27 -#ifndef PREFLOW_PUSH_MAX_FLOW_H
18.28 -#define PREFLOW_PUSH_MAX_FLOW_H
18.29 -
18.30 -#define A 1
18.31 -
18.32 -#include <algorithm>
18.33 -#include <vector>
18.34 -#include <stack>
18.35 -
18.36 -#include <reverse_bfs.h>
18.37 -
18.38 -
18.39 -namespace hugo {
18.40 -
18.41 - template <typename Graph, typename T,
18.42 - typename FlowMap=typename Graph::EdgeMap<T>, typename CapMap=typename Graph::EdgeMap<T>,
18.43 - typename IntMap=typename Graph::NodeMap<int>, typename TMap=typename Graph::NodeMap<T> >
18.44 - class preflow_push_max_flow {
18.45 -
18.46 - typedef typename Graph::NodeIt NodeIt;
18.47 - typedef typename Graph::EachNodeIt EachNodeIt;
18.48 - typedef typename Graph::OutEdgeIt OutEdgeIt;
18.49 - typedef typename Graph::InEdgeIt InEdgeIt;
18.50 -
18.51 - Graph& G;
18.52 - NodeIt s;
18.53 - NodeIt t;
18.54 - IntMap level;
18.55 - CapMap& capacity;
18.56 - int empty_level; //an empty level in the end of run()
18.57 - T value;
18.58 -
18.59 - public:
18.60 -
18.61 - preflow_push_max_flow(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
18.62 - G(_G), s(_s), t(_t), level(_G), capacity(_capacity) { }
18.63 -
18.64 -
18.65 - /*
18.66 - The run() function runs a modified version of the
18.67 - highest label preflow-push, which only
18.68 - finds a maximum preflow, hence giving the value of a maximum flow.
18.69 - */
18.70 - void run() {
18.71 -
18.72 - int n=G.nodeNum();
18.73 - int b=n-2;
18.74 - /*
18.75 - b is a bound on the highest level of an active node.
18.76 - */
18.77 -
18.78 - IntMap level(G,n);
18.79 - TMap excess(G);
18.80 - FlowMap flow(G,0);
18.81 -
18.82 - std::vector<int> numb(n);
18.83 - /*
18.84 - The number of nodes on level i < n. It is
18.85 - initialized to n+1, because of the reverse_bfs-part.
18.86 - */
18.87 -
18.88 - std::vector<std::stack<NodeIt> > stack(n);
18.89 - //Stack of the active nodes in level i.
18.90 -
18.91 -
18.92 - /*Reverse_bfs from t, to find the starting level.*/
18.93 - level.set(t,0);
18.94 - std::queue<NodeIt> bfs_queue;
18.95 - bfs_queue.push(t);
18.96 -
18.97 - while (!bfs_queue.empty()) {
18.98 -
18.99 - NodeIt v=bfs_queue.front();
18.100 - bfs_queue.pop();
18.101 - int l=level.get(v)+1;
18.102 -
18.103 - for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
18.104 - NodeIt w=G.tail(e);
18.105 - if ( level.get(w) == n ) {
18.106 - bfs_queue.push(w);
18.107 - ++numb[l];
18.108 - level.set(w, l);
18.109 - }
18.110 - }
18.111 - }
18.112 -
18.113 - level.set(s,n);
18.114 -
18.115 -
18.116 -
18.117 - /* Starting flow. It is everywhere 0 at the moment. */
18.118 - for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
18.119 - {
18.120 - if ( capacity.get(e) == 0 ) continue;
18.121 - NodeIt w=G.head(e);
18.122 - if ( level.get(w) < n ) {
18.123 - if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w);
18.124 - flow.set(e, capacity.get(e));
18.125 - excess.set(w, excess.get(w)+capacity.get(e));
18.126 - }
18.127 - }
18.128 -
18.129 - /*
18.130 - End of preprocessing
18.131 - */
18.132 -
18.133 -
18.134 - /*
18.135 - Push/relabel on the highest level active nodes.
18.136 - */
18.137 - /*While there exists an active node.*/
18.138 - while (b) {
18.139 - if ( stack[b].empty() ) {
18.140 - --b;
18.141 - continue;
18.142 - }
18.143 -
18.144 - NodeIt w=stack[b].top(); //w is a highest label active node.
18.145 - stack[b].pop();
18.146 - int lev=level.get(w);
18.147 - int exc=excess.get(w);
18.148 - int newlevel=2*n-2; //In newlevel we bound the next level of w.
18.149 -
18.150 - // if ( level.get(w) < n ) { //Nem tudom ez mukodik-e
18.151 - for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
18.152 -
18.153 - if ( flow.get(e) == capacity.get(e) ) continue;
18.154 - NodeIt v=G.head(e);
18.155 - //e=wv
18.156 -
18.157 - if( lev > level.get(v) ) {
18.158 - /*Push is allowed now*/
18.159 -
18.160 - if ( excess.get(v)==0 && v != s && v !=t )
18.161 - stack[level.get(v)].push(v);
18.162 - /*v becomes active.*/
18.163 -
18.164 - int cap=capacity.get(e);
18.165 - int flo=flow.get(e);
18.166 - int remcap=cap-flo;
18.167 -
18.168 - if ( remcap >= exc ) {
18.169 - /*A nonsaturating push.*/
18.170 -
18.171 - flow.set(e, flo+exc);
18.172 - excess.set(v, excess.get(v)+exc);
18.173 - exc=0;
18.174 - break;
18.175 -
18.176 - } else {
18.177 - /*A saturating push.*/
18.178 -
18.179 - flow.set(e, cap );
18.180 - excess.set(v, excess.get(v)+remcap);
18.181 - exc-=remcap;
18.182 - }
18.183 - } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
18.184 -
18.185 - } //for out edges wv
18.186 -
18.187 -
18.188 - if ( exc > 0 ) {
18.189 - for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
18.190 -
18.191 - if( flow.get(e) == 0 ) continue;
18.192 - NodeIt v=G.tail(e);
18.193 - //e=vw
18.194 -
18.195 - if( lev > level.get(v) ) {
18.196 - /*Push is allowed now*/
18.197 -
18.198 - if ( excess.get(v)==0 && v != s && v !=t)
18.199 - stack[level.get(v)].push(v);
18.200 - /*v becomes active.*/
18.201 -
18.202 - int flo=flow.get(e);
18.203 -
18.204 - if ( flo >= exc ) {
18.205 - /*A nonsaturating push.*/
18.206 -
18.207 - flow.set(e, flo-exc);
18.208 - excess.set(v, excess.get(v)+exc);
18.209 - exc=0;
18.210 - break;
18.211 - } else {
18.212 - /*A saturating push.*/
18.213 -
18.214 - excess.set(v, excess.get(v)+flo);
18.215 - exc-=flo;
18.216 - flow.set(e,0);
18.217 - }
18.218 - } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
18.219 -
18.220 - } //for in edges vw
18.221 -
18.222 - } // if w still has excess after the out edge for cycle
18.223 -
18.224 - excess.set(w, exc);
18.225 -
18.226 -
18.227 - /*
18.228 - Relabel
18.229 - */
18.230 -
18.231 - if ( exc > 0 ) {
18.232 - //now 'lev' is the old level of w
18.233 - level.set(w,++newlevel);
18.234 - --numb[lev];
18.235 -
18.236 - if ( !numb[lev] && lev < A*n ) { //If the level of w gets empty.
18.237 -
18.238 - for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
18.239 - if (level.get(v) > lev ) level.set(v,n);
18.240 - }
18.241 - for (int i=lev+1 ; i!=n ; ++i) numb[i]=0;
18.242 - if ( newlevel < n ) newlevel=n;
18.243 - } else if ( newlevel < n ) {
18.244 - ++numb[newlevel];
18.245 - stack[newlevel].push(w);
18.246 - b=newlevel;
18.247 - }
18.248 - }
18.249 -
18.250 -
18.251 -
18.252 - } //while(b)
18.253 -
18.254 - value=excess.get(t);
18.255 - /*Max flow value.*/
18.256 -
18.257 -
18.258 - /*
18.259 - We count empty_level. The nodes above this level is a mincut.
18.260 - */
18.261 - while(true) {
18.262 - if(numb[empty_level]) ++empty_level;
18.263 - else break;
18.264 - }
18.265 -
18.266 - } // void run()
18.267 -
18.268 -
18.269 -
18.270 - /*
18.271 - Returns the maximum value of a flow.
18.272 - */
18.273 -
18.274 - T maxflow() {
18.275 - return value;
18.276 - }
18.277 -
18.278 -
18.279 -
18.280 - /*
18.281 - Returns a minimum cut.
18.282 - */
18.283 - template<typename CutMap>
18.284 - void mincut(CutMap& M) {
18.285 -
18.286 - for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v) {
18.287 - if ( level.get(v) > empty_level ) M.set(v, true);
18.288 - }
18.289 - }
18.290 -
18.291 -
18.292 - };
18.293 -}//namespace hugo
18.294 -#endif
18.295 -
18.296 -
18.297 -
18.298 -
18.299 -
19.1 --- a/src/work/jacint/preflow_push_max_flow.hh Fri Feb 20 22:01:02 2004 +0000
19.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
19.3 @@ -1,315 +0,0 @@
19.4 -/*
19.5 -preflow_push_max_flow_hh
19.6 -by jacint.
19.7 -Runs a preflow push algorithm with the modification,
19.8 -that we do not push on nodes with level at least n.
19.9 -Moreover, if a level gets empty, we put all nodes above that
19.10 -level to level n. Hence, in the end, we arrive at a maximum preflow
19.11 -with value of a max flow value. An empty level gives a minimum cut.
19.12 -
19.13 -Member functions:
19.14 -
19.15 -void run() : runs the algorithm
19.16 -
19.17 - The following functions should be used after run() was already run.
19.18 -
19.19 -T maxflow() : returns the value of a maximum flow
19.20 -
19.21 -node_property_vector<graph_type, bool> mincut(): returns a
19.22 - characteristic vector of a minimum cut.
19.23 -*/
19.24 -
19.25 -#ifndef PREFLOW_PUSH_MAX_FLOW_HH
19.26 -#define PREFLOW_PUSH_MAX_FLOW_HH
19.27 -
19.28 -#include <algorithm>
19.29 -#include <vector>
19.30 -#include <stack>
19.31 -
19.32 -#include <marci_list_graph.hh>
19.33 -#include <marci_graph_traits.hh>
19.34 -#include <marci_property_vector.hh>
19.35 -#include <reverse_bfs.hh>
19.36 -
19.37 -
19.38 -namespace hugo {
19.39 -
19.40 - template <typename graph_type, typename T>
19.41 - class preflow_push_max_flow {
19.42 -
19.43 - typedef typename graph_traits<graph_type>::node_iterator node_iterator;
19.44 - typedef typename graph_traits<graph_type>::each_node_iterator each_node_iterator;
19.45 - typedef typename graph_traits<graph_type>::out_edge_iterator out_edge_iterator;
19.46 - typedef typename graph_traits<graph_type>::in_edge_iterator in_edge_iterator;
19.47 -
19.48 - graph_type& G;
19.49 - node_iterator s;
19.50 - node_iterator t;
19.51 - edge_property_vector<graph_type, T>& capacity;
19.52 - T value;
19.53 - node_property_vector<graph_type, bool> mincutvector;
19.54 -
19.55 -
19.56 -
19.57 - public:
19.58 -
19.59 - preflow_push_max_flow(graph_type& _G, node_iterator _s, node_iterator _t, edge_property_vector<graph_type, T>& _capacity) : G(_G), s(_s), t(_t), capacity(_capacity), mincutvector(_G, false) { }
19.60 -
19.61 -
19.62 - /*
19.63 - The run() function runs a modified version of the highest label preflow-push, which only
19.64 - finds a maximum preflow, hence giving the value of a maximum flow.
19.65 - */
19.66 - void run() {
19.67 -
19.68 - edge_property_vector<graph_type, T> flow(G, 0); //the flow value, 0 everywhere
19.69 - node_property_vector<graph_type, int> level(G); //level of node
19.70 - node_property_vector<graph_type, T> excess(G); //excess of node
19.71 -
19.72 - int n=number_of(G.first_node()); //number of nodes
19.73 - int b=n-2;
19.74 - /*b is a bound on the highest level of an active node. In the beginning it is at most n-2.*/
19.75 -
19.76 - std::vector<int> numb(n); //The number of nodes on level i < n.
19.77 -
19.78 - std::vector<std::stack<node_iterator> > stack(2*n-1); //Stack of the active nodes in level i.
19.79 -
19.80 -
19.81 -
19.82 - /*Reverse_bfs from t, to find the starting level.*/
19.83 -
19.84 - reverse_bfs<list_graph> bfs(G, t);
19.85 - bfs.run();
19.86 - for(each_node_iterator v=G.first_node(); v.valid(); ++v)
19.87 - {
19.88 - int dist=bfs.dist(v);
19.89 - level.put(v, dist);
19.90 - ++numb[dist];
19.91 - }
19.92 -
19.93 - /*The level of s is fixed to n*/
19.94 - level.put(s,n);
19.95 -
19.96 -
19.97 - /* Starting flow. It is everywhere 0 at the moment. */
19.98 -
19.99 - for(out_edge_iterator i=G.first_out_edge(s); i.valid(); ++i)
19.100 - {
19.101 - node_iterator w=G.head(i);
19.102 - flow.put(i, capacity.get(i));
19.103 - stack[bfs.dist(w)].push(w);
19.104 - excess.put(w, capacity.get(i));
19.105 - }
19.106 -
19.107 -
19.108 - /*
19.109 - End of preprocessing
19.110 - */
19.111 -
19.112 -
19.113 -
19.114 -
19.115 - /*
19.116 - Push/relabel on the highest level active nodes.
19.117 - */
19.118 -
19.119 - /*While there exists an active node.*/
19.120 - while (b) {
19.121 -
19.122 - /*We decrease the bound if there is no active node of level b.*/
19.123 - if (stack[b].empty()) {
19.124 - --b;
19.125 - } else {
19.126 -
19.127 - node_iterator w=stack[b].top(); //w is the highest label active node.
19.128 - stack[b].pop(); //We delete w from the stack.
19.129 -
19.130 - int newlevel=2*n-2; //In newlevel we maintain the next level of w.
19.131 -
19.132 - for(out_edge_iterator e=G.first_out_edge(w); e.valid(); ++e) {
19.133 - node_iterator v=G.head(e);
19.134 - /*e is the edge wv.*/
19.135 -
19.136 - if (flow.get(e)<capacity.get(e)) {
19.137 - /*e is an edge of the residual graph */
19.138 -
19.139 - if(level.get(w)==level.get(v)+1) {
19.140 - /*Push is allowed now*/
19.141 -
19.142 - if (capacity.get(e)-flow.get(e) > excess.get(w)) {
19.143 - /*A nonsaturating push.*/
19.144 -
19.145 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
19.146 - /*v becomes active.*/
19.147 -
19.148 - flow.put(e, flow.get(e)+excess.get(w));
19.149 - excess.put(v, excess.get(v)+excess.get(w));
19.150 - excess.put(w,0);
19.151 - //std::cout << w << " " << v <<" elore elen nonsat pump " << std::endl;
19.152 - break;
19.153 - } else {
19.154 - /*A saturating push.*/
19.155 -
19.156 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
19.157 - /*v becomes active.*/
19.158 -
19.159 - excess.put(v, excess.get(v)+capacity.get(e)-flow.get(e));
19.160 - excess.put(w, excess.get(w)-capacity.get(e)+flow.get(e));
19.161 - flow.put(e, capacity.get(e));
19.162 - //std::cout << w <<" " << v <<" elore elen sat pump " << std::endl;
19.163 - if (excess.get(w)==0) break;
19.164 - /*If w is not active any more, then we go on to the next node.*/
19.165 -
19.166 - } // if (capacity.get(e)-flow.get(e) > excess.get(w))
19.167 - } // if (level.get(w)==level.get(v)+1)
19.168 -
19.169 - else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);}
19.170 -
19.171 - } //if (flow.get(e)<capacity.get(e))
19.172 -
19.173 - } //for(out_edge_iterator e=G.first_out_edge(w); e.valid(); ++e)
19.174 -
19.175 -
19.176 -
19.177 - for(in_edge_iterator e=G.first_in_edge(w); e.valid(); ++e) {
19.178 - node_iterator v=G.tail(e);
19.179 - /*e is the edge vw.*/
19.180 -
19.181 - if (excess.get(w)==0) break;
19.182 - /*It may happen, that w became inactive in the first 'for' cycle.*/
19.183 -
19.184 - if(flow.get(e)>0) {
19.185 - /*e is an edge of the residual graph */
19.186 -
19.187 - if(level.get(w)==level.get(v)+1) {
19.188 - /*Push is allowed now*/
19.189 -
19.190 - if (flow.get(e) > excess.get(w)) {
19.191 - /*A nonsaturating push.*/
19.192 -
19.193 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
19.194 - /*v becomes active.*/
19.195 -
19.196 - flow.put(e, flow.get(e)-excess.get(w));
19.197 - excess.put(v, excess.get(v)+excess.get(w));
19.198 - excess.put(w,0);
19.199 - //std::cout << v << " " << w << " vissza elen nonsat pump " << std::endl;
19.200 - break;
19.201 - } else {
19.202 - /*A saturating push.*/
19.203 -
19.204 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
19.205 - /*v becomes active.*/
19.206 -
19.207 - flow.put(e,0);
19.208 - excess.put(v, excess.get(v)+flow.get(e));
19.209 - excess.put(w, excess.get(w)-flow.get(e));
19.210 - //std::cout << v <<" " << w << " vissza elen sat pump " << std::endl;
19.211 - if (excess.get(w)==0) { break;}
19.212 - } //if (flow.get(e) > excess.get(v))
19.213 - } //if(level.get(w)==level.get(v)+1)
19.214 -
19.215 - else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);}
19.216 - //std::cout << "Leveldecrease of node " << w << " to " << newlevel << std::endl;
19.217 -
19.218 - } //if (flow.get(e)>0)
19.219 -
19.220 - } //for in-edge
19.221 -
19.222 -
19.223 -
19.224 -
19.225 - /*
19.226 - Relabel
19.227 - */
19.228 - if (excess.get(w)>0) {
19.229 - /*Now newlevel <= n*/
19.230 -
19.231 - int l=level.get(w); //l is the old level of w.
19.232 - --numb[l];
19.233 -
19.234 - if (newlevel == n) {
19.235 - level.put(w,n);
19.236 -
19.237 - } else {
19.238 -
19.239 - if (numb[l]) {
19.240 - /*If the level of w remains nonempty.*/
19.241 -
19.242 - level.put(w,++newlevel);
19.243 - ++numb[newlevel];
19.244 - stack[newlevel].push(w);
19.245 - b=newlevel;
19.246 - } else {
19.247 - /*If the level of w gets empty.*/
19.248 -
19.249 - for (each_node_iterator v=G.first_node() ; v.valid() ; ++v) {
19.250 - if (level.get(v) >= l ) {
19.251 - level.put(v,n);
19.252 - }
19.253 - }
19.254 -
19.255 - for (int i=l+1 ; i!=n ; ++i) numb[i]=0;
19.256 - } //if (numb[l])
19.257 -
19.258 - } // if (newlevel = n)
19.259 -
19.260 - } // if (excess.get(w)>0)
19.261 -
19.262 -
19.263 - } //else
19.264 -
19.265 - } //while(b)
19.266 -
19.267 - value=excess.get(t);
19.268 - /*Max flow value.*/
19.269 -
19.270 -
19.271 -
19.272 - /*
19.273 - We find an empty level, e. The nodes above this level give
19.274 - a minimum cut.
19.275 - */
19.276 -
19.277 - int e=1;
19.278 -
19.279 - while(e) {
19.280 - if(numb[e]) ++e;
19.281 - else break;
19.282 - }
19.283 - for (each_node_iterator v=G.first_node(); v.valid(); ++v) {
19.284 - if (level.get(v) > e) mincutvector.put(v, true);
19.285 - }
19.286 -
19.287 -
19.288 - } // void run()
19.289 -
19.290 -
19.291 -
19.292 - /*
19.293 - Returns the maximum value of a flow.
19.294 - */
19.295 -
19.296 - T maxflow() {
19.297 - return value;
19.298 - }
19.299 -
19.300 -
19.301 -
19.302 - /*
19.303 - Returns a minimum cut.
19.304 - */
19.305 -
19.306 - node_property_vector<graph_type, bool> mincut() {
19.307 - return mincutvector;
19.308 - }
19.309 -
19.310 -
19.311 - };
19.312 -}//namespace hugo
19.313 -#endif
19.314 -
19.315 -
19.316 -
19.317 -
19.318 -