1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/jacint/preflow.cc Tue Apr 27 22:59:15 2004 +0000
1.3 @@ -0,0 +1,239 @@
1.4 +#include <iostream>
1.5 +
1.6 +#include <smart_graph.h>
1.7 +#include <dimacs.h>
1.8 +#include <preflow.h>
1.9 +#include <time_measure.h>
1.10 +
1.11 +using namespace hugo;
1.12 +
1.13 +int main(int, char **) {
1.14 +
1.15 + typedef SmartGraph Graph;
1.16 +
1.17 + typedef Graph::Node Node;
1.18 + typedef Graph::EdgeIt EdgeIt;
1.19 +
1.20 + Graph G;
1.21 + Node s, t;
1.22 + Graph::EdgeMap<int> cap(G);
1.23 + readDimacsMaxFlow(std::cin, G, s, t, cap);
1.24 + Timer ts;
1.25 +
1.26 + std::cout <<
1.27 + "\n Testing preflow.h on a graph with " <<
1.28 + G.nodeNum() << " nodes and " << G.edgeNum() << " edges..."
1.29 + << std::endl;
1.30 +
1.31 +
1.32 + Graph::EdgeMap<int> flow(G,0);
1.33 + Preflow<Graph, int> preflow_test(G, s, t, cap, flow);
1.34 + std::cout << "\nCalling run() (flow must be constant zero)..."<<std::endl;
1.35 + ts.reset();
1.36 + preflow_test.run();
1.37 + std::cout << "Elapsed time: " << ts << std::endl;
1.38 +
1.39 + Graph::NodeMap<bool> mincut(G);
1.40 + preflow_test.minMinCut(mincut);
1.41 + int min_min_cut_value=0;
1.42 + Graph::NodeMap<bool> cut(G);
1.43 + preflow_test.minCut(cut);
1.44 + int min_cut_value=0;
1.45 + Graph::NodeMap<bool> maxcut(G);
1.46 + preflow_test.maxMinCut(maxcut);
1.47 + int max_min_cut_value=0;
1.48 + EdgeIt e;
1.49 + for(G.first(e); G.valid(e); G.next(e)) {
1.50 + int c=cap[e];
1.51 + if (mincut[G.tail(e)] && !mincut[G.head(e)]) min_min_cut_value+=c;
1.52 + if (cut[G.tail(e)] && !cut[G.head(e)]) min_cut_value+=c;
1.53 + if (maxcut[G.tail(e)] && !maxcut[G.head(e)]) max_min_cut_value+=c;
1.54 + }
1.55 +
1.56 + std::cout << "\nChecking the result: " <<std::endl;
1.57 + std::cout << "Flow value: "<< preflow_test.flowValue() << std::endl;
1.58 + std::cout << "Min cut value: "<< min_cut_value << std::endl;
1.59 + std::cout << "Min min cut value: "<< min_min_cut_value << std::endl;
1.60 + std::cout << "Max min cut value: "<< max_min_cut_value <<
1.61 + std::endl;
1.62 +
1.63 + if ( preflow_test.flowValue() == min_cut_value &&
1.64 + min_cut_value == min_min_cut_value &&
1.65 + min_min_cut_value == max_min_cut_value )
1.66 + std::cout << "They are equal. " <<std::endl;
1.67 +
1.68 +
1.69 +
1.70 +
1.71 +
1.72 + Preflow<Graph, int> preflow_test2(G, s, t, cap, flow);
1.73 + std::cout << "\n\nCalling preflow(GEN_FLOW) with the given maximum flow..."<<std::endl;
1.74 + ts.reset();
1.75 + preflow_test2.preflow(preflow_test2.GEN_FLOW);
1.76 + std::cout << "Elapsed time: " << ts << std::endl;
1.77 +
1.78 + Graph::NodeMap<bool> mincut2(G);
1.79 + preflow_test.minMinCut(mincut2);
1.80 + int min_min_cut2_value=0;
1.81 + Graph::NodeMap<bool> cut2(G);
1.82 + preflow_test.minCut(cut2);
1.83 + int min_cut2_value=0;
1.84 + Graph::NodeMap<bool> maxcut2(G);
1.85 + preflow_test.maxMinCut(maxcut2);
1.86 + int max_min_cut2_value=0;
1.87 + for(G.first(e); G.valid(e); G.next(e)) {
1.88 + int c=cap[e];
1.89 + if (mincut2[G.tail(e)] && !mincut2[G.head(e)]) min_min_cut2_value+=c;
1.90 + if (cut2[G.tail(e)] && !cut2[G.head(e)]) min_cut2_value+=c;
1.91 + if (maxcut2[G.tail(e)] && !maxcut2[G.head(e)]) max_min_cut2_value+=c;
1.92 + }
1.93 +
1.94 + std::cout << "\nThe given flow value is "
1.95 + << preflow_test2.flowValue();
1.96 +
1.97 + if ( preflow_test2.flowValue() == min_cut2_value &&
1.98 + min_cut2_value == min_min_cut2_value &&
1.99 + min_min_cut2_value == max_min_cut2_value )
1.100 + std::cout <<", which is equal to all three min cut values."
1.101 + <<std::endl;
1.102 +
1.103 +
1.104 +
1.105 +
1.106 +
1.107 + Graph::EdgeMap<int> flow3(G,0);
1.108 + Preflow<Graph, int> preflow_test3(G, s, t, cap, flow3);
1.109 + std::cout << "\n\nCalling preflowPhase0(PREFLOW) on the constant zero flow..."<<std::endl;
1.110 + ts.reset();
1.111 + preflow_test3.preflowPhase0(preflow_test3.PREFLOW);
1.112 + std::cout << "Elapsed time: " << ts << std::endl;
1.113 + Graph::NodeMap<bool> actcut3(G);
1.114 + std::cout << "\nCalling actMinCut()..."<<std::endl;
1.115 + preflow_test3.actMinCut(actcut3);
1.116 + std::cout << "Calling preflowPhase1() on the given flow..."<<std::endl;
1.117 + ts.reset();
1.118 + preflow_test3.preflowPhase1();
1.119 + std::cout << "Elapsed time: " << ts << std::endl;
1.120 +
1.121 + int act_min_cut3_value=0;
1.122 +
1.123 + Graph::NodeMap<bool> mincut3(G);
1.124 + preflow_test.minMinCut(mincut3);
1.125 + int min_min_cut3_value=0;
1.126 +
1.127 + Graph::NodeMap<bool> cut3(G);
1.128 + preflow_test.minCut(cut3);
1.129 + int min_cut3_value=0;
1.130 +
1.131 + Graph::NodeMap<bool> maxcut3(G);
1.132 + preflow_test.maxMinCut(maxcut3);
1.133 + int max_min_cut3_value=0;
1.134 +
1.135 + for(G.first(e); G.valid(e); G.next(e)) {
1.136 + int c=cap[e];
1.137 + if (mincut3[G.tail(e)] && !mincut3[G.head(e)]) min_min_cut3_value+=c;
1.138 + if (cut3[G.tail(e)] && !cut3[G.head(e)]) min_cut3_value+=c;
1.139 + if (maxcut3[G.tail(e)] && !maxcut3[G.head(e)]) max_min_cut3_value+=c;
1.140 + if (actcut3[G.tail(e)] && !actcut3[G.head(e)]) act_min_cut3_value+=c;
1.141 + }
1.142 +
1.143 + std::cout << "\nThe min cut value given by actMinCut() after phase 0 is "<<
1.144 + act_min_cut3_value;
1.145 +
1.146 + if ( preflow_test3.flowValue() == min_cut3_value &&
1.147 + min_cut3_value == min_min_cut3_value &&
1.148 + min_min_cut3_value == max_min_cut3_value &&
1.149 + max_min_cut3_value == act_min_cut3_value ) {
1.150 + std::cout <<
1.151 + ", which is equal to the given flow value and to all three min cut values after phase 1."
1.152 + <<std::endl;
1.153 + }
1.154 +
1.155 +
1.156 +
1.157 +
1.158 +
1.159 + Graph::EdgeMap<int> flow4(G,0);
1.160 + Preflow<Graph, int> preflow_test4(G, s, t, cap, flow4);
1.161 + std::cout <<
1.162 + "\n\nCalling preflow(PREFLOW) with the constant 0 flow, the result is f..."
1.163 + <<std::endl;
1.164 + preflow_test4.preflow(preflow_test4.PREFLOW);
1.165 +
1.166 + std::cout << "Swapping the source and the target, "<<std::endl;
1.167 + std::cout << "by calling resetSource(t) and resetTarget(s)..."
1.168 + <<std::endl;
1.169 + preflow_test4.resetSource(t);
1.170 + preflow_test4.resetTarget(s);
1.171 +
1.172 + std::cout <<
1.173 + "Calling preflow(PREFLOW) to find a maximum t-s flow starting with flow f..."
1.174 + <<std::endl;
1.175 + preflow_test4.preflow(preflow_test4.PREFLOW);
1.176 +
1.177 + Graph::NodeMap<bool> mincut4(G);
1.178 + preflow_test4.minMinCut(mincut4);
1.179 + int min_min_cut4_value=0;
1.180 + Graph::NodeMap<bool> cut4(G);
1.181 + preflow_test4.minCut(cut4);
1.182 + int min_cut4_value=0;
1.183 + Graph::NodeMap<bool> maxcut4(G);
1.184 + preflow_test4.maxMinCut(maxcut4);
1.185 + int max_min_cut4_value=0;
1.186 + for(G.first(e); G.valid(e); G.next(e)) {
1.187 + int c=cap[e];
1.188 + if (mincut4[G.tail(e)] && !mincut4[G.head(e)]) min_min_cut4_value+=c;
1.189 + if (cut4[G.tail(e)] && !cut4[G.head(e)]) min_cut4_value+=c;
1.190 + if (maxcut4[G.tail(e)] && !maxcut4[G.head(e)]) max_min_cut4_value+=c;
1.191 + }
1.192 +
1.193 + std::cout << "\nThe given flow value is "
1.194 + << preflow_test4.flowValue();
1.195 +
1.196 + if ( preflow_test4.flowValue() == min_cut4_value &&
1.197 + min_cut4_value == min_min_cut4_value &&
1.198 + min_min_cut4_value == max_min_cut4_value )
1.199 + std::cout <<", which is equal to all three min cut values."
1.200 + <<std::endl;
1.201 +
1.202 +
1.203 +
1.204 +
1.205 + Graph::EdgeMap<int> flow5(G,0);
1.206 + std::cout << "Resetting the stored flow to constant zero, by calling resetFlow..."
1.207 + <<std::endl;
1.208 + preflow_test4.resetFlow(flow5);
1.209 + std::cout <<
1.210 + "Calling preflow(GEN_FLOW) to find a maximum t-s flow "<<std::endl;
1.211 + std::cout <<
1.212 + "starting with this constant zero flow..." <<std::endl;
1.213 + preflow_test4.preflow(preflow_test4.GEN_FLOW);
1.214 +
1.215 + Graph::NodeMap<bool> mincut5(G);
1.216 + preflow_test4.minMinCut(mincut5);
1.217 + int min_min_cut5_value=0;
1.218 + Graph::NodeMap<bool> cut5(G);
1.219 + preflow_test4.minCut(cut5);
1.220 + int min_cut5_value=0;
1.221 + Graph::NodeMap<bool> maxcut5(G);
1.222 + preflow_test4.maxMinCut(maxcut5);
1.223 + int max_min_cut5_value=0;
1.224 + for(G.first(e); G.valid(e); G.next(e)) {
1.225 + int c=cap[e];
1.226 + if (mincut5[G.tail(e)] && !mincut5[G.head(e)]) min_min_cut5_value+=c;
1.227 + if (cut5[G.tail(e)] && !cut5[G.head(e)]) min_cut5_value+=c;
1.228 + if (maxcut5[G.tail(e)] && !maxcut5[G.head(e)]) max_min_cut5_value+=c;
1.229 + }
1.230 +
1.231 + std::cout << "\nThe given flow value is "
1.232 + << preflow_test4.flowValue();
1.233 +
1.234 + if ( preflow_test4.flowValue() == min_cut5_value &&
1.235 + min_cut5_value == min_min_cut5_value &&
1.236 + min_min_cut5_value == max_min_cut5_value )
1.237 + std::cout <<", which is equal to all three min cut values."
1.238 + <<std::endl<<std::endl;
1.239 +
1.240 +
1.241 + return 0;
1.242 +}
2.1 --- a/src/work/jacint/preflow.h Tue Apr 27 22:29:11 2004 +0000
2.2 +++ b/src/work/jacint/preflow.h Tue Apr 27 22:59:15 2004 +0000
2.3 @@ -1,20 +1,15 @@
2.4 // -*- C++ -*-
2.5
2.6 -//run gyorsan tudna adni a minmincutot a 2 fazis elejen , ne vegyuk be konstruktorba egy cutmapet?
2.7 -//constzero jo igy?
2.8 -
2.9 -//majd marci megmondja betegyem-e bfs-t meg resgraphot
2.10 -
2.11 /*
2.12 Heuristics:
2.13 2 phase
2.14 gap
2.15 list 'level_list' on the nodes on level i implemented by hand
2.16 - stack 'active' on the active nodes on level i implemented by hand
2.17 + stack 'active' on the active nodes on level i
2.18 runs heuristic 'highest label' for H1*n relabels
2.19 runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
2.20
2.21 -Parameters H0 and H1 are initialized to 20 and 10.
2.22 +Parameters H0 and H1 are initialized to 20 and 1.
2.23
2.24 Constructors:
2.25
2.26 @@ -28,7 +23,7 @@
2.27 T flowValue() : returns the value of a maximum flow
2.28
2.29 void minMinCut(CutMap& M) : sets M to the characteristic vector of the
2.30 - minimum min cut. M should be a map of bools initialized to false.
2.31 + minimum min cut. M should be a map of bools initialized to false. ??Is it OK?
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 @@ -36,8 +31,6 @@
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 -FIXME reset
2.40 -
2.41 */
2.42
2.43 #ifndef HUGO_PREFLOW_H
2.44 @@ -48,6 +41,7 @@
2.45
2.46 #include <vector>
2.47 #include <queue>
2.48 +#include <stack>
2.49
2.50 namespace hugo {
2.51
2.52 @@ -57,494 +51,223 @@
2.53 class Preflow {
2.54
2.55 typedef typename Graph::Node Node;
2.56 - typedef typename Graph::Edge Edge;
2.57 typedef typename Graph::NodeIt NodeIt;
2.58 typedef typename Graph::OutEdgeIt OutEdgeIt;
2.59 typedef typename Graph::InEdgeIt InEdgeIt;
2.60 -
2.61 +
2.62 + typedef typename std::vector<std::stack<Node> > VecStack;
2.63 + typedef typename Graph::template NodeMap<Node> NNMap;
2.64 + typedef typename std::vector<Node> VecNode;
2.65 +
2.66 const Graph& G;
2.67 Node s;
2.68 Node t;
2.69 - const CapMap& capacity;
2.70 - FlowMap& flow;
2.71 - T value;
2.72 - bool constzero;
2.73 + CapMap* capacity;
2.74 + FlowMap* flow;
2.75 + int n; //the number of nodes of G
2.76 + typename Graph::template NodeMap<int> level;
2.77 + typename Graph::template NodeMap<T> excess;
2.78 +
2.79
2.80 public:
2.81 +
2.82 + enum flowEnum{
2.83 + ZERO_FLOW=0,
2.84 + GEN_FLOW=1,
2.85 + PREFLOW=2
2.86 + };
2.87 +
2.88 Preflow(Graph& _G, Node _s, Node _t, CapMap& _capacity,
2.89 - FlowMap& _flow, bool _constzero ) :
2.90 - G(_G), s(_s), t(_t), capacity(_capacity), flow(_flow), constzero(_constzero) {}
2.91 + FlowMap& _flow) :
2.92 + G(_G), s(_s), t(_t), capacity(&_capacity),
2.93 + flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
2.94 +
2.95 + void run() {
2.96 + preflow( ZERO_FLOW );
2.97 + }
2.98
2.99 -
2.100 - void run() {
2.101 + void preflow( flowEnum fe ) {
2.102 + preflowPhase0(fe);
2.103 + preflowPhase1();
2.104 + }
2.105 +
2.106 + void preflowPhase0( flowEnum fe ) {
2.107
2.108 - value=0; //for the subsequent runs
2.109 -
2.110 - bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd
2.111 - int n=G.nodeNum();
2.112 int heur0=(int)(H0*n); //time while running 'bound decrease'
2.113 int heur1=(int)(H1*n); //time while running 'highest label'
2.114 int heur=heur1; //starting time interval (#of relabels)
2.115 + int numrelabel=0;
2.116 +
2.117 bool what_heur=1;
2.118 - /*
2.119 - what_heur is 0 in case 'bound decrease'
2.120 - and 1 in case 'highest label'
2.121 - */
2.122 + //It is 0 in case 'bound decrease' and 1 in case 'highest label'
2.123 +
2.124 bool end=false;
2.125 - /*
2.126 - Needed for 'bound decrease', 'true'
2.127 - means no active nodes are above bound b.
2.128 - */
2.129 - int relabel=0;
2.130 + //Needed for 'bound decrease', true means no active nodes are above bound b.
2.131 +
2.132 int k=n-2; //bound on the highest level under n containing a node
2.133 int b=k; //bound on the highest level under n of an active node
2.134
2.135 - typename Graph::template NodeMap<int> level(G,n);
2.136 - typename Graph::template NodeMap<T> excess(G);
2.137 + VecStack active(n);
2.138 +
2.139 + NNMap left(G,INVALID);
2.140 + NNMap right(G,INVALID);
2.141 + VecNode level_list(n,INVALID);
2.142 + //List of the nodes in level i<n, set to n.
2.143
2.144 - std::vector<Node> active(n-1,INVALID);
2.145 - typename Graph::template NodeMap<Node> next(G,INVALID);
2.146 - //Stack of the active nodes in level i < n.
2.147 - //We use it in both phases.
2.148 -
2.149 - typename Graph::template NodeMap<Node> left(G,INVALID);
2.150 - typename Graph::template NodeMap<Node> right(G,INVALID);
2.151 - std::vector<Node> level_list(n,INVALID);
2.152 - /*
2.153 - List of the nodes in level i<n.
2.154 - */
2.155 -
2.156 -
2.157 - if ( constzero ) {
2.158 -
2.159 - /*Reverse_bfs from t, to find the starting level.*/
2.160 - level.set(t,0);
2.161 - std::queue<Node> bfs_queue;
2.162 - bfs_queue.push(t);
2.163 -
2.164 - while (!bfs_queue.empty()) {
2.165 + NodeIt v;
2.166 + for(G.first(v); G.valid(v); G.next(v)) level.set(v,n);
2.167 + //setting each node to level n
2.168 +
2.169 + switch ( fe ) {
2.170 + case PREFLOW:
2.171 + {
2.172 + //counting the excess
2.173 + NodeIt v;
2.174 + for(G.first(v); G.valid(v); G.next(v)) {
2.175 + T exc=0;
2.176
2.177 - Node v=bfs_queue.front();
2.178 - bfs_queue.pop();
2.179 - int l=level[v]+1;
2.180 + InEdgeIt e;
2.181 + for(G.first(e,v); G.valid(e); G.next(e)) exc+=(*flow)[e];
2.182 + OutEdgeIt f;
2.183 + for(G.first(f,v); G.valid(f); G.next(f)) exc-=(*flow)[f];
2.184 +
2.185 + excess.set(v,exc);
2.186 +
2.187 + //putting the active nodes into the stack
2.188 + int lev=level[v];
2.189 + if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
2.190 + }
2.191 + break;
2.192 + }
2.193 + case GEN_FLOW:
2.194 + {
2.195 + //Counting the excess of t
2.196 + T exc=0;
2.197
2.198 InEdgeIt e;
2.199 - for(G.first(e,v); G.valid(e); G.next(e)) {
2.200 - Node w=G.tail(e);
2.201 - if ( level[w] == n && w != s ) {
2.202 - bfs_queue.push(w);
2.203 - Node first=level_list[l];
2.204 - if ( G.valid(first) ) left.set(first,w);
2.205 - right.set(w,first);
2.206 - level_list[l]=w;
2.207 - level.set(w, l);
2.208 + for(G.first(e,t); G.valid(e); G.next(e)) exc+=(*flow)[e];
2.209 + OutEdgeIt f;
2.210 + for(G.first(f,t); G.valid(f); G.next(f)) exc-=(*flow)[f];
2.211 +
2.212 + excess.set(t,exc);
2.213 +
2.214 + break;
2.215 + }
2.216 + }
2.217 +
2.218 + preflowPreproc( fe, active, level_list, left, right );
2.219 + //End of preprocessing
2.220 +
2.221 +
2.222 + //Push/relabel on the highest level active nodes.
2.223 + while ( true ) {
2.224 + if ( b == 0 ) {
2.225 + if ( !what_heur && !end && k > 0 ) {
2.226 + b=k;
2.227 + end=true;
2.228 + } else break;
2.229 + }
2.230 +
2.231 + if ( active[b].empty() ) --b;
2.232 + else {
2.233 + end=false;
2.234 + Node w=active[b].top();
2.235 + active[b].pop();
2.236 + int newlevel=push(w,active);
2.237 + if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list,
2.238 + left, right, b, k, what_heur);
2.239 +
2.240 + ++numrelabel;
2.241 + if ( numrelabel >= heur ) {
2.242 + numrelabel=0;
2.243 + if ( what_heur ) {
2.244 + what_heur=0;
2.245 + heur=heur0;
2.246 + end=false;
2.247 + } else {
2.248 + what_heur=1;
2.249 + heur=heur1;
2.250 + b=k;
2.251 }
2.252 }
2.253 + }
2.254 + }
2.255 + }
2.256 +
2.257 +
2.258 + void preflowPhase1() {
2.259 +
2.260 + int k=n-2; //bound on the highest level under n containing a node
2.261 + int b=k; //bound on the highest level under n of an active node
2.262 +
2.263 + VecStack active(n);
2.264 + level.set(s,0);
2.265 + std::queue<Node> bfs_queue;
2.266 + bfs_queue.push(s);
2.267 +
2.268 + while (!bfs_queue.empty()) {
2.269 +
2.270 + Node v=bfs_queue.front();
2.271 + bfs_queue.pop();
2.272 + int l=level[v]+1;
2.273 +
2.274 + InEdgeIt e;
2.275 + for(G.first(e,v); G.valid(e); G.next(e)) {
2.276 + if ( (*capacity)[e] == (*flow)[e] ) continue;
2.277 + Node u=G.tail(e);
2.278 + if ( level[u] >= n ) {
2.279 + bfs_queue.push(u);
2.280 + level.set(u, l);
2.281 + if ( excess[u] > 0 ) active[l].push(u);
2.282 + }
2.283 }
2.284 -
2.285 - //the starting flow
2.286 - OutEdgeIt e;
2.287 - for(G.first(e,s); G.valid(e); G.next(e))
2.288 - {
2.289 - T c=capacity[e];
2.290 - if ( c == 0 ) continue;
2.291 - Node w=G.head(e);
2.292 - if ( level[w] < n ) {
2.293 - if ( excess[w] == 0 && w!=t ) {
2.294 - next.set(w,active[level[w]]);
2.295 - active[level[w]]=w;
2.296 - }
2.297 - flow.set(e, c);
2.298 - excess.set(w, excess[w]+c);
2.299 +
2.300 + OutEdgeIt f;
2.301 + for(G.first(f,v); G.valid(f); G.next(f)) {
2.302 + if ( 0 == (*flow)[f] ) continue;
2.303 + Node u=G.head(f);
2.304 + if ( level[u] >= n ) {
2.305 + bfs_queue.push(u);
2.306 + level.set(u, l);
2.307 + if ( excess[u] > 0 ) active[l].push(u);
2.308 }
2.309 }
2.310 }
2.311 - else
2.312 - {
2.313 -
2.314 - /*
2.315 - Reverse_bfs from t in the residual graph,
2.316 - to find the starting level.
2.317 - */
2.318 - level.set(t,0);
2.319 - std::queue<Node> bfs_queue;
2.320 - bfs_queue.push(t);
2.321 -
2.322 - while (!bfs_queue.empty()) {
2.323 -
2.324 - Node v=bfs_queue.front();
2.325 - bfs_queue.pop();
2.326 - int l=level[v]+1;
2.327 -
2.328 - InEdgeIt e;
2.329 - for(G.first(e,v); G.valid(e); G.next(e)) {
2.330 - if ( capacity[e] == flow[e] ) continue;
2.331 - Node w=G.tail(e);
2.332 - if ( level[w] == n && w != s ) {
2.333 - bfs_queue.push(w);
2.334 - Node first=level_list[l];
2.335 - if ( G.valid(first) ) left.set(first,w);
2.336 - right.set(w,first);
2.337 - level_list[l]=w;
2.338 - level.set(w, l);
2.339 - }
2.340 - }
2.341 -
2.342 - OutEdgeIt f;
2.343 - for(G.first(f,v); G.valid(f); G.next(f)) {
2.344 - if ( 0 == flow[f] ) continue;
2.345 - Node w=G.head(f);
2.346 - if ( level[w] == n && w != s ) {
2.347 - bfs_queue.push(w);
2.348 - Node first=level_list[l];
2.349 - if ( G.valid(first) ) left.set(first,w);
2.350 - right.set(w,first);
2.351 - level_list[l]=w;
2.352 - level.set(w, l);
2.353 - }
2.354 - }
2.355 - }
2.356 -
2.357 -
2.358 - /*
2.359 - Counting the excess
2.360 - */
2.361 - NodeIt v;
2.362 - for(G.first(v); G.valid(v); G.next(v)) {
2.363 - T exc=0;
2.364 + b=n-2;
2.365
2.366 - InEdgeIt e;
2.367 - for(G.first(e,v); G.valid(e); G.next(e)) exc+=flow[e];
2.368 - OutEdgeIt f;
2.369 - for(G.first(f,v); G.valid(f); G.next(f)) exc-=flow[e];
2.370 -
2.371 - excess.set(v,exc);
2.372 -
2.373 - //putting the active nodes into the stack
2.374 - int lev=level[v];
2.375 - if ( exc > 0 && lev < n ) {
2.376 - next.set(v,active[lev]);
2.377 - active[lev]=v;
2.378 - }
2.379 - }
2.380 -
2.381 -
2.382 - //the starting flow
2.383 - OutEdgeIt e;
2.384 - for(G.first(e,s); G.valid(e); G.next(e))
2.385 - {
2.386 - T rem=capacity[e]-flow[e];
2.387 - if ( rem == 0 ) continue;
2.388 - Node w=G.head(e);
2.389 - if ( level[w] < n ) {
2.390 - if ( excess[w] == 0 && w!=t ) {
2.391 - next.set(w,active[level[w]]);
2.392 - active[level[w]]=w;
2.393 - }
2.394 - flow.set(e, capacity[e]);
2.395 - excess.set(w, excess[w]+rem);
2.396 - }
2.397 - }
2.398 -
2.399 - InEdgeIt f;
2.400 - for(G.first(f,s); G.valid(f); G.next(f))
2.401 - {
2.402 - if ( flow[f] == 0 ) continue;
2.403 - Node w=G.head(f);
2.404 - if ( level[w] < n ) {
2.405 - if ( excess[w] == 0 && w!=t ) {
2.406 - next.set(w,active[level[w]]);
2.407 - active[level[w]]=w;
2.408 - }
2.409 - excess.set(w, excess[w]+flow[f]);
2.410 - flow.set(f, 0);
2.411 - }
2.412 - }
2.413 - }
2.414 -
2.415 -
2.416 -
2.417 -
2.418 - /*
2.419 - End of preprocessing
2.420 - */
2.421 -
2.422 -
2.423 -
2.424 - /*
2.425 - Push/relabel on the highest level active nodes.
2.426 - */
2.427 while ( true ) {
2.428
2.429 - if ( b == 0 ) {
2.430 - if ( phase ) break;
2.431 -
2.432 - if ( !what_heur && !end && k > 0 ) {
2.433 - b=k;
2.434 - end=true;
2.435 - } else {
2.436 - phase=1;
2.437 - level.set(s,0);
2.438 - std::queue<Node> bfs_queue;
2.439 - bfs_queue.push(s);
2.440 -
2.441 - while (!bfs_queue.empty()) {
2.442 -
2.443 - Node v=bfs_queue.front();
2.444 - bfs_queue.pop();
2.445 - int l=level[v]+1;
2.446 -
2.447 - InEdgeIt e;
2.448 - for(G.first(e,v); G.valid(e); G.next(e)) {
2.449 - if ( capacity[e] == flow[e] ) continue;
2.450 - Node u=G.tail(e);
2.451 - if ( level[u] >= n ) {
2.452 - bfs_queue.push(u);
2.453 - level.set(u, l);
2.454 - if ( excess[u] > 0 ) {
2.455 - next.set(u,active[l]);
2.456 - active[l]=u;
2.457 - }
2.458 - }
2.459 - }
2.460 -
2.461 - OutEdgeIt f;
2.462 - for(G.first(f,v); G.valid(f); G.next(f)) {
2.463 - if ( 0 == flow[f] ) continue;
2.464 - Node u=G.head(f);
2.465 - if ( level[u] >= n ) {
2.466 - bfs_queue.push(u);
2.467 - level.set(u, l);
2.468 - if ( excess[u] > 0 ) {
2.469 - next.set(u,active[l]);
2.470 - active[l]=u;
2.471 - }
2.472 - }
2.473 - }
2.474 - }
2.475 - b=n-2;
2.476 - }
2.477 -
2.478 - }
2.479 -
2.480 -
2.481 - if ( !G.valid(active[b]) ) --b;
2.482 + if ( b == 0 ) break;
2.483 +
2.484 + if ( active[b].empty() ) --b;
2.485 else {
2.486 - end=false;
2.487 + Node w=active[b].top();
2.488 + active[b].pop();
2.489 + int newlevel=push(w,active);
2.490
2.491 - Node w=active[b];
2.492 - active[b]=next[w];
2.493 - int lev=level[w];
2.494 - T exc=excess[w];
2.495 - int newlevel=n; //bound on the next level of w
2.496 -
2.497 - OutEdgeIt e;
2.498 - for(G.first(e,w); G.valid(e); G.next(e)) {
2.499 -
2.500 - if ( flow[e] == capacity[e] ) continue;
2.501 - Node v=G.head(e);
2.502 - //e=wv
2.503 -
2.504 - if( lev > level[v] ) {
2.505 - /*Push is allowed now*/
2.506 -
2.507 - if ( excess[v]==0 && v!=t && v!=s ) {
2.508 - int lev_v=level[v];
2.509 - next.set(v,active[lev_v]);
2.510 - active[lev_v]=v;
2.511 - }
2.512 -
2.513 - T cap=capacity[e];
2.514 - T flo=flow[e];
2.515 - T remcap=cap-flo;
2.516 -
2.517 - if ( remcap >= exc ) {
2.518 - /*A nonsaturating push.*/
2.519 -
2.520 - flow.set(e, flo+exc);
2.521 - excess.set(v, excess[v]+exc);
2.522 - exc=0;
2.523 - break;
2.524 -
2.525 - } else {
2.526 - /*A saturating push.*/
2.527 -
2.528 - flow.set(e, cap);
2.529 - excess.set(v, excess[v]+remcap);
2.530 - exc-=remcap;
2.531 - }
2.532 - } else if ( newlevel > level[v] ){
2.533 - newlevel = level[v];
2.534 - }
2.535 -
2.536 - } //for out edges wv
2.537 -
2.538 -
2.539 - if ( exc > 0 ) {
2.540 - InEdgeIt e;
2.541 - for(G.first(e,w); G.valid(e); G.next(e)) {
2.542 -
2.543 - if( flow[e] == 0 ) continue;
2.544 - Node v=G.tail(e);
2.545 - //e=vw
2.546 -
2.547 - if( lev > level[v] ) {
2.548 - /*Push is allowed now*/
2.549 -
2.550 - if ( excess[v]==0 && v!=t && v!=s ) {
2.551 - int lev_v=level[v];
2.552 - next.set(v,active[lev_v]);
2.553 - active[lev_v]=v;
2.554 - }
2.555 -
2.556 - T flo=flow[e];
2.557 -
2.558 - if ( flo >= exc ) {
2.559 - /*A nonsaturating push.*/
2.560 -
2.561 - flow.set(e, flo-exc);
2.562 - excess.set(v, excess[v]+exc);
2.563 - exc=0;
2.564 - break;
2.565 - } else {
2.566 - /*A saturating push.*/
2.567 -
2.568 - excess.set(v, excess[v]+flo);
2.569 - exc-=flo;
2.570 - flow.set(e,0);
2.571 - }
2.572 - } else if ( newlevel > level[v] ) {
2.573 - newlevel = level[v];
2.574 - }
2.575 - } //for in edges vw
2.576 -
2.577 - } // if w still has excess after the out edge for cycle
2.578 -
2.579 - excess.set(w, exc);
2.580 -
2.581 - /*
2.582 - Relabel
2.583 - */
2.584 -
2.585 -
2.586 - if ( exc > 0 ) {
2.587 - //now 'lev' is the old level of w
2.588 -
2.589 - if ( phase ) {
2.590 + //relabel
2.591 + if ( excess[w] > 0 ) {
2.592 level.set(w,++newlevel);
2.593 - next.set(w,active[newlevel]);
2.594 - active[newlevel]=w;
2.595 + active[newlevel].push(w);
2.596 b=newlevel;
2.597 - } else {
2.598 - //unlacing starts
2.599 - Node right_n=right[w];
2.600 - Node left_n=left[w];
2.601 -
2.602 - if ( G.valid(right_n) ) {
2.603 - if ( G.valid(left_n) ) {
2.604 - right.set(left_n, right_n);
2.605 - left.set(right_n, left_n);
2.606 - } else {
2.607 - level_list[lev]=right_n;
2.608 - left.set(right_n, INVALID);
2.609 - }
2.610 - } else {
2.611 - if ( G.valid(left_n) ) {
2.612 - right.set(left_n, INVALID);
2.613 - } else {
2.614 - level_list[lev]=INVALID;
2.615 - }
2.616 - }
2.617 - //unlacing ends
2.618 -
2.619 - if ( !G.valid(level_list[lev]) ) {
2.620 -
2.621 - //gapping starts
2.622 - for (int i=lev; i!=k ; ) {
2.623 - Node v=level_list[++i];
2.624 - while ( G.valid(v) ) {
2.625 - level.set(v,n);
2.626 - v=right[v];
2.627 - }
2.628 - level_list[i]=INVALID;
2.629 - if ( !what_heur ) active[i]=INVALID;
2.630 - }
2.631 -
2.632 - level.set(w,n);
2.633 - b=lev-1;
2.634 - k=b;
2.635 - //gapping ends
2.636 -
2.637 - } else {
2.638 -
2.639 - if ( newlevel == n ) level.set(w,n);
2.640 - else {
2.641 - level.set(w,++newlevel);
2.642 - next.set(w,active[newlevel]);
2.643 - active[newlevel]=w;
2.644 - if ( what_heur ) b=newlevel;
2.645 - if ( k < newlevel ) ++k; //now k=newlevel
2.646 - Node first=level_list[newlevel];
2.647 - if ( G.valid(first) ) left.set(first,w);
2.648 - right.set(w,first);
2.649 - left.set(w,INVALID);
2.650 - level_list[newlevel]=w;
2.651 - }
2.652 - }
2.653 -
2.654 -
2.655 - ++relabel;
2.656 - if ( relabel >= heur ) {
2.657 - relabel=0;
2.658 - if ( what_heur ) {
2.659 - what_heur=0;
2.660 - heur=heur0;
2.661 - end=false;
2.662 - } else {
2.663 - what_heur=1;
2.664 - heur=heur1;
2.665 - b=k;
2.666 - }
2.667 - }
2.668 - } //phase 0
2.669 -
2.670 -
2.671 - } // if ( exc > 0 )
2.672 -
2.673 -
2.674 + }
2.675 } // if stack[b] is nonempty
2.676 -
2.677 } // while(true)
2.678 -
2.679 -
2.680 - value = excess[t];
2.681 - /*Max flow value.*/
2.682 -
2.683 - } //void run()
2.684 -
2.685 -
2.686 -
2.687 -
2.688 -
2.689 - /*
2.690 - Returns the maximum value of a flow.
2.691 - */
2.692 -
2.693 - T flowValue() {
2.694 - return value;
2.695 }
2.696
2.697
2.698 - FlowMap Flow() {
2.699 - return flow;
2.700 - }
2.701 + //Returns the maximum value of a flow.
2.702 + T flowValue() {
2.703 + return excess[t];
2.704 + }
2.705
2.706 -
2.707 -
2.708 - void Flow(FlowMap& _flow ) {
2.709 + //should be used only between preflowPhase0 and preflowPhase1
2.710 + template<typename _CutMap>
2.711 + void actMinCut(_CutMap& M) {
2.712 NodeIt v;
2.713 - for(G.first(v) ; G.valid(v); G.next(v))
2.714 - _flow.set(v,flow[v]);
2.715 + for(G.first(v); G.valid(v); G.next(v))
2.716 + if ( level[v] < n ) M.set(v,false);
2.717 + else M.set(v,true);
2.718 }
2.719
2.720
2.721 @@ -552,7 +275,6 @@
2.722 /*
2.723 Returns the minimum min cut, by a bfs from s in the residual graph.
2.724 */
2.725 -
2.726 template<typename _CutMap>
2.727 void minMinCut(_CutMap& M) {
2.728
2.729 @@ -568,7 +290,7 @@
2.730 OutEdgeIt e;
2.731 for(G.first(e,w) ; G.valid(e); G.next(e)) {
2.732 Node v=G.head(e);
2.733 - if (!M[v] && flow[e] < capacity[e] ) {
2.734 + if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
2.735 queue.push(v);
2.736 M.set(v, true);
2.737 }
2.738 @@ -577,7 +299,7 @@
2.739 InEdgeIt f;
2.740 for(G.first(f,w) ; G.valid(f); G.next(f)) {
2.741 Node v=G.tail(f);
2.742 - if (!M[v] && flow[f] > 0 ) {
2.743 + if (!M[v] && (*flow)[f] > 0 ) {
2.744 queue.push(v);
2.745 M.set(v, true);
2.746 }
2.747 @@ -594,10 +316,15 @@
2.748
2.749 template<typename _CutMap>
2.750 void maxMinCut(_CutMap& M) {
2.751 -
2.752 +
2.753 + NodeIt v;
2.754 + for(G.first(v) ; G.valid(v); G.next(v)) {
2.755 + M.set(v, true);
2.756 + }
2.757 +
2.758 std::queue<Node> queue;
2.759
2.760 - M.set(t,true);
2.761 + M.set(t,false);
2.762 queue.push(t);
2.763
2.764 while (!queue.empty()) {
2.765 @@ -608,56 +335,319 @@
2.766 InEdgeIt e;
2.767 for(G.first(e,w) ; G.valid(e); G.next(e)) {
2.768 Node v=G.tail(e);
2.769 - if (!M[v] && flow[e] < capacity[e] ) {
2.770 + if (M[v] && (*flow)[e] < (*capacity)[e] ) {
2.771 queue.push(v);
2.772 - M.set(v, true);
2.773 + M.set(v, false);
2.774 }
2.775 }
2.776
2.777 OutEdgeIt f;
2.778 for(G.first(f,w) ; G.valid(f); G.next(f)) {
2.779 Node v=G.head(f);
2.780 - if (!M[v] && flow[f] > 0 ) {
2.781 + if (M[v] && (*flow)[f] > 0 ) {
2.782 queue.push(v);
2.783 - M.set(v, true);
2.784 + M.set(v, false);
2.785 }
2.786 }
2.787 }
2.788 -
2.789 - NodeIt v;
2.790 - for(G.first(v) ; G.valid(v); G.next(v)) {
2.791 - M.set(v, !M[v]);
2.792 - }
2.793 -
2.794 }
2.795
2.796
2.797 -
2.798 template<typename CutMap>
2.799 void minCut(CutMap& M) {
2.800 minMinCut(M);
2.801 }
2.802
2.803
2.804 - void reset_target (Node _t) {t=_t;}
2.805 - void reset_source (Node _s) {s=_s;}
2.806 + void resetTarget (const Node _t) {t=_t;}
2.807 +
2.808 + void resetSource (const Node _s) {s=_s;}
2.809
2.810 - template<typename _CapMap>
2.811 - void reset_cap (_CapMap _cap) {capacity=_cap;}
2.812 -
2.813 - template<typename _FlowMap>
2.814 - void reset_cap (_FlowMap _flow, bool _constzero) {
2.815 - flow=_flow;
2.816 - constzero=_constzero;
2.817 + void resetCap (const CapMap& _cap) {
2.818 + capacity=&_cap;
2.819 + }
2.820 +
2.821 + void resetFlow (FlowMap& _flow) {
2.822 + flow=&_flow;
2.823 }
2.824
2.825
2.826 + private:
2.827 +
2.828 + int push(const Node w, VecStack& active) {
2.829 +
2.830 + int lev=level[w];
2.831 + T exc=excess[w];
2.832 + int newlevel=n; //bound on the next level of w
2.833 +
2.834 + OutEdgeIt e;
2.835 + for(G.first(e,w); G.valid(e); G.next(e)) {
2.836 +
2.837 + if ( (*flow)[e] == (*capacity)[e] ) continue;
2.838 + Node v=G.head(e);
2.839 +
2.840 + if( lev > level[v] ) { //Push is allowed now
2.841 +
2.842 + if ( excess[v]==0 && v!=t && v!=s ) {
2.843 + int lev_v=level[v];
2.844 + active[lev_v].push(v);
2.845 + }
2.846 +
2.847 + T cap=(*capacity)[e];
2.848 + T flo=(*flow)[e];
2.849 + T remcap=cap-flo;
2.850 +
2.851 + if ( remcap >= exc ) { //A nonsaturating push.
2.852 +
2.853 + flow->set(e, flo+exc);
2.854 + excess.set(v, excess[v]+exc);
2.855 + exc=0;
2.856 + break;
2.857 +
2.858 + } else { //A saturating push.
2.859 + flow->set(e, cap);
2.860 + excess.set(v, excess[v]+remcap);
2.861 + exc-=remcap;
2.862 + }
2.863 + } else if ( newlevel > level[v] ) newlevel = level[v];
2.864 + } //for out edges wv
2.865 +
2.866 + if ( exc > 0 ) {
2.867 + InEdgeIt e;
2.868 + for(G.first(e,w); G.valid(e); G.next(e)) {
2.869 +
2.870 + if( (*flow)[e] == 0 ) continue;
2.871 + Node v=G.tail(e);
2.872 +
2.873 + if( lev > level[v] ) { //Push is allowed now
2.874 +
2.875 + if ( excess[v]==0 && v!=t && v!=s ) {
2.876 + int lev_v=level[v];
2.877 + active[lev_v].push(v);
2.878 + }
2.879 +
2.880 + T flo=(*flow)[e];
2.881 +
2.882 + if ( flo >= exc ) { //A nonsaturating push.
2.883 +
2.884 + flow->set(e, flo-exc);
2.885 + excess.set(v, excess[v]+exc);
2.886 + exc=0;
2.887 + break;
2.888 + } else { //A saturating push.
2.889 +
2.890 + excess.set(v, excess[v]+flo);
2.891 + exc-=flo;
2.892 + flow->set(e,0);
2.893 + }
2.894 + } else if ( newlevel > level[v] ) newlevel = level[v];
2.895 + } //for in edges vw
2.896 +
2.897 + } // if w still has excess after the out edge for cycle
2.898 +
2.899 + excess.set(w, exc);
2.900 +
2.901 + return newlevel;
2.902 + }
2.903 +
2.904 +
2.905 + void preflowPreproc ( flowEnum fe, VecStack& active,
2.906 + VecNode& level_list, NNMap& left, NNMap& right ) {
2.907 +
2.908 + std::queue<Node> bfs_queue;
2.909 +
2.910 + switch ( fe ) {
2.911 + case ZERO_FLOW:
2.912 + {
2.913 + //Reverse_bfs from t, to find the starting level.
2.914 + level.set(t,0);
2.915 + bfs_queue.push(t);
2.916 +
2.917 + while (!bfs_queue.empty()) {
2.918 +
2.919 + Node v=bfs_queue.front();
2.920 + bfs_queue.pop();
2.921 + int l=level[v]+1;
2.922 +
2.923 + InEdgeIt e;
2.924 + for(G.first(e,v); G.valid(e); G.next(e)) {
2.925 + Node w=G.tail(e);
2.926 + if ( level[w] == n && w != s ) {
2.927 + bfs_queue.push(w);
2.928 + Node first=level_list[l];
2.929 + if ( G.valid(first) ) left.set(first,w);
2.930 + right.set(w,first);
2.931 + level_list[l]=w;
2.932 + level.set(w, l);
2.933 + }
2.934 + }
2.935 + }
2.936 +
2.937 + //the starting flow
2.938 + OutEdgeIt e;
2.939 + for(G.first(e,s); G.valid(e); G.next(e))
2.940 + {
2.941 + T c=(*capacity)[e];
2.942 + if ( c == 0 ) continue;
2.943 + Node w=G.head(e);
2.944 + if ( level[w] < n ) {
2.945 + if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
2.946 + flow->set(e, c);
2.947 + excess.set(w, excess[w]+c);
2.948 + }
2.949 + }
2.950 + break;
2.951 + }
2.952 +
2.953 + case GEN_FLOW:
2.954 + case PREFLOW:
2.955 + {
2.956 + //Reverse_bfs from t in the residual graph,
2.957 + //to find the starting level.
2.958 + level.set(t,0);
2.959 + bfs_queue.push(t);
2.960 +
2.961 + while (!bfs_queue.empty()) {
2.962 +
2.963 + Node v=bfs_queue.front();
2.964 + bfs_queue.pop();
2.965 + int l=level[v]+1;
2.966 +
2.967 + InEdgeIt e;
2.968 + for(G.first(e,v); G.valid(e); G.next(e)) {
2.969 + if ( (*capacity)[e] == (*flow)[e] ) continue;
2.970 + Node w=G.tail(e);
2.971 + if ( level[w] == n && w != s ) {
2.972 + bfs_queue.push(w);
2.973 + Node first=level_list[l];
2.974 + if ( G.valid(first) ) left.set(first,w);
2.975 + right.set(w,first);
2.976 + level_list[l]=w;
2.977 + level.set(w, l);
2.978 + }
2.979 + }
2.980 +
2.981 + OutEdgeIt f;
2.982 + for(G.first(f,v); G.valid(f); G.next(f)) {
2.983 + if ( 0 == (*flow)[f] ) continue;
2.984 + Node w=G.head(f);
2.985 + if ( level[w] == n && w != s ) {
2.986 + bfs_queue.push(w);
2.987 + Node first=level_list[l];
2.988 + if ( G.valid(first) ) left.set(first,w);
2.989 + right.set(w,first);
2.990 + level_list[l]=w;
2.991 + level.set(w, l);
2.992 + }
2.993 + }
2.994 + }
2.995 +
2.996 +
2.997 + //the starting flow
2.998 + OutEdgeIt e;
2.999 + for(G.first(e,s); G.valid(e); G.next(e))
2.1000 + {
2.1001 + T rem=(*capacity)[e]-(*flow)[e];
2.1002 + if ( rem == 0 ) continue;
2.1003 + Node w=G.head(e);
2.1004 + if ( level[w] < n ) {
2.1005 + if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
2.1006 + flow->set(e, (*capacity)[e]);
2.1007 + excess.set(w, excess[w]+rem);
2.1008 + }
2.1009 + }
2.1010 +
2.1011 + InEdgeIt f;
2.1012 + for(G.first(f,s); G.valid(f); G.next(f))
2.1013 + {
2.1014 + if ( (*flow)[f] == 0 ) continue;
2.1015 + Node w=G.tail(f);
2.1016 + if ( level[w] < n ) {
2.1017 + if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
2.1018 + excess.set(w, excess[w]+(*flow)[f]);
2.1019 + flow->set(f, 0);
2.1020 + }
2.1021 + }
2.1022 + break;
2.1023 + } //case PREFLOW
2.1024 + }
2.1025 + } //preflowPreproc
2.1026 +
2.1027 +
2.1028 +
2.1029 + void relabel( const Node w, int newlevel, VecStack& active,
2.1030 + VecNode& level_list, NNMap& left,
2.1031 + NNMap& right, int& b, int& k, const bool what_heur ) {
2.1032 +
2.1033 + T lev=level[w];
2.1034 +
2.1035 + Node right_n=right[w];
2.1036 + Node left_n=left[w];
2.1037 +
2.1038 + //unlacing starts
2.1039 + if ( G.valid(right_n) ) {
2.1040 + if ( G.valid(left_n) ) {
2.1041 + right.set(left_n, right_n);
2.1042 + left.set(right_n, left_n);
2.1043 + } else {
2.1044 + level_list[lev]=right_n;
2.1045 + left.set(right_n, INVALID);
2.1046 + }
2.1047 + } else {
2.1048 + if ( G.valid(left_n) ) {
2.1049 + right.set(left_n, INVALID);
2.1050 + } else {
2.1051 + level_list[lev]=INVALID;
2.1052 + }
2.1053 + }
2.1054 + //unlacing ends
2.1055 +
2.1056 + if ( !G.valid(level_list[lev]) ) {
2.1057 +
2.1058 + //gapping starts
2.1059 + for (int i=lev; i!=k ; ) {
2.1060 + Node v=level_list[++i];
2.1061 + while ( G.valid(v) ) {
2.1062 + level.set(v,n);
2.1063 + v=right[v];
2.1064 + }
2.1065 + level_list[i]=INVALID;
2.1066 + if ( !what_heur ) {
2.1067 + while ( !active[i].empty() ) {
2.1068 + active[i].pop(); //FIXME: ezt szebben kene
2.1069 + }
2.1070 + }
2.1071 + }
2.1072 +
2.1073 + level.set(w,n);
2.1074 + b=lev-1;
2.1075 + k=b;
2.1076 + //gapping ends
2.1077 +
2.1078 + } else {
2.1079 +
2.1080 + if ( newlevel == n ) level.set(w,n);
2.1081 + else {
2.1082 + level.set(w,++newlevel);
2.1083 + active[newlevel].push(w);
2.1084 + if ( what_heur ) b=newlevel;
2.1085 + if ( k < newlevel ) ++k; //now k=newlevel
2.1086 + Node first=level_list[newlevel];
2.1087 + if ( G.valid(first) ) left.set(first,w);
2.1088 + right.set(w,first);
2.1089 + left.set(w,INVALID);
2.1090 + level_list[newlevel]=w;
2.1091 + }
2.1092 + }
2.1093 +
2.1094 + } //relabel
2.1095 +
2.1096
2.1097 };
2.1098
2.1099 } //namespace hugo
2.1100
2.1101 -#endif //PREFLOW_H
2.1102 +#endif //HUGO_PREFLOW_H
2.1103
2.1104
2.1105
3.1 --- a/src/work/jacint/preflow_res_comp.cc Tue Apr 27 22:29:11 2004 +0000
3.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
3.3 @@ -1,125 +0,0 @@
3.4 -/*
3.5 -The only difference between preflow.h and preflow_res.h is that the latter
3.6 -uses the ResGraphWrapper, while the first does not. (Bfs is implemented by
3.7 -hand in both.) This test program runs Preflow and PreflowRes on the same
3.8 -graph, tests the result of these implementations and writes the running time
3.9 -of them. */
3.10 -#include <iostream>
3.11 -
3.12 -#include <smart_graph.h>
3.13 -#include <dimacs.h>
3.14 -#include <preflow.h>
3.15 -#include <preflow_res.h>
3.16 -#include <time_measure.h>
3.17 -
3.18 -using namespace hugo;
3.19 -
3.20 -int main(int, char **) {
3.21 -
3.22 - typedef SmartGraph Graph;
3.23 -
3.24 - typedef Graph::Node Node;
3.25 - typedef Graph::EdgeIt EdgeIt;
3.26 -
3.27 - Graph G;
3.28 - Node s, t;
3.29 - Graph::EdgeMap<int> cap(G);
3.30 - readDimacsMaxFlow(std::cin, G, s, t, cap);
3.31 - Timer ts;
3.32 -
3.33 - std::cout <<
3.34 - "\n In which way are we faster: using ResGraphWrapper or not?"
3.35 - <<std::endl;
3.36 - std::cout <<
3.37 - "\n Running preflow.h on a graph with " <<
3.38 - G.nodeNum() << " nodes and " << G.edgeNum() << " edges..."
3.39 - << std::endl<<std::endl;
3.40 -
3.41 -
3.42 -
3.43 - Graph::EdgeMap<int> flow(G,0);
3.44 - Preflow<Graph, int> max_flow_test(G, s, t, cap, flow, 1);
3.45 - ts.reset();
3.46 - max_flow_test.run();
3.47 - std::cout << "Elapsed time NOT using the ResGraphWrapper: " << std::endl
3.48 - <<ts << std::endl;
3.49 -
3.50 - Graph::NodeMap<bool> mincut(G);
3.51 - max_flow_test.minMinCut(mincut);
3.52 - int min_min_cut_value=0;
3.53 - EdgeIt e;
3.54 - for(G.first(e); G.valid(e); G.next(e)) {
3.55 - if (mincut[G.tail(e)] && !mincut[G.head(e)]) min_min_cut_value+=cap[e];
3.56 - }
3.57 -
3.58 - Graph::NodeMap<bool> cut(G);
3.59 - max_flow_test.minCut(cut);
3.60 - int min_cut_value=0;
3.61 - for(G.first(e); G.valid(e); G.next(e)) {
3.62 - if (cut[G.tail(e)] && !cut[G.head(e)])
3.63 - min_cut_value+=cap[e];
3.64 - }
3.65 -
3.66 - Graph::NodeMap<bool> maxcut(G);
3.67 - max_flow_test.maxMinCut(maxcut);
3.68 - int max_min_cut_value=0;
3.69 - for(G.first(e); G.valid(e); G.next(e)) {
3.70 - if (maxcut[G.tail(e)] && !maxcut[G.head(e)])
3.71 - max_min_cut_value+=cap[e];
3.72 - }
3.73 -
3.74 - std::cout << "\n Checking the result: " <<std::endl;
3.75 - std::cout << "Flow value: "<< max_flow_test.flowValue() << std::endl;
3.76 - std::cout << "Min cut value: "<< min_cut_value << std::endl;
3.77 - std::cout << "Min min cut value: "<< min_min_cut_value << std::endl;
3.78 - std::cout << "Max min cut value: "<< max_min_cut_value <<
3.79 - std::endl;
3.80 -
3.81 - if ( max_flow_test.flowValue() == min_cut_value &&
3.82 - min_cut_value == min_min_cut_value &&
3.83 - min_min_cut_value == max_min_cut_value )
3.84 - std::cout << "They are equal! " <<std::endl<< std::endl<<"\n";
3.85 -
3.86 - Graph::EdgeMap<int> flow2(G,0);
3.87 - PreflowRes<Graph, int> max_flow_test2(G, s, t, cap, flow2, 1);
3.88 - ts.reset();
3.89 - max_flow_test2.run();
3.90 - std::cout << "Elapsed time using the ResGraphWrapper: " << std::endl
3.91 - << ts << std::endl;
3.92 -
3.93 - Graph::NodeMap<bool> mincut2(G);
3.94 - max_flow_test2.minMinCut(mincut2);
3.95 - int min_min_cut_value2=0;
3.96 - for(G.first(e); G.valid(e); G.next(e)) {
3.97 - if (mincut2[G.tail(e)] && !mincut2[G.head(e)]) min_min_cut_value2+=cap[e];
3.98 - }
3.99 -
3.100 - Graph::NodeMap<bool> cut2(G);
3.101 - max_flow_test2.minCut(cut2);
3.102 - int min_cut_value2=0;
3.103 - for(G.first(e); G.valid(e); G.next(e)) {
3.104 - if (cut2[G.tail(e)] && !cut2[G.head(e)])
3.105 - min_cut_value2+=cap[e];
3.106 - }
3.107 -
3.108 - Graph::NodeMap<bool> maxcut2(G);
3.109 - max_flow_test2.maxMinCut(maxcut2);
3.110 - int max_min_cut_value2=0;
3.111 - for(G.first(e); G.valid(e); G.next(e)) {
3.112 - if (maxcut2[G.tail(e)] && !maxcut2[G.head(e)])
3.113 - max_min_cut_value2+=cap[e];
3.114 - }
3.115 -
3.116 - std::cout << "\n Checking the result: " <<std::endl;
3.117 - std::cout << "Flow value: "<< max_flow_test2.flowValue() << std::endl;
3.118 - std::cout << "Min cut value: "<< min_cut_value2 << std::endl;
3.119 - std::cout << "Min min cut value: "<< min_min_cut_value2 << std::endl;
3.120 - std::cout << "Max min cut value: "<< max_min_cut_value2 <<
3.121 - std::endl;
3.122 - if ( max_flow_test.flowValue() == min_cut_value &&
3.123 - min_cut_value == min_min_cut_value &&
3.124 - min_min_cut_value == max_min_cut_value )
3.125 - std::cout << "They are equal! " <<std::endl<<"/n";
3.126 -
3.127 - return 0;
3.128 -}