1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/jacint/preflow_push_max_flow.hh Fri Jan 30 14:53:17 2004 +0000
1.3 @@ -0,0 +1,315 @@
1.4 +/*
1.5 +preflow_push_max_flow_hh
1.6 +by jacint.
1.7 +Runs a preflow push algorithm with the modification,
1.8 +that we do not push on nodes with level at least n.
1.9 +Moreover, if a level gets empty, we put all nodes above that
1.10 +level to level n. Hence, in the end, we arrive at a maximum preflow
1.11 +with value of a max flow value. An empty level gives a minimum cut.
1.12 +
1.13 +Member functions:
1.14 +
1.15 +void run() : runs the algorithm
1.16 +
1.17 + The following functions should be used after run() was already run.
1.18 +
1.19 +T maxflow() : returns the value of a maximum flow
1.20 +
1.21 +node_property_vector<graph_type, bool> mincut(): returns a
1.22 + characteristic vector of a minimum cut.
1.23 +*/
1.24 +
1.25 +#ifndef PREFLOW_PUSH_MAX_FLOW_HH
1.26 +#define PREFLOW_PUSH_MAX_FLOW_HH
1.27 +
1.28 +#include <algorithm>
1.29 +#include <vector>
1.30 +#include <stack>
1.31 +
1.32 +#include <marci_list_graph.hh>
1.33 +#include <marci_graph_traits.hh>
1.34 +#include <marci_property_vector.hh>
1.35 +#include <reverse_bfs.hh>
1.36 +
1.37 +
1.38 +namespace marci {
1.39 +
1.40 + template <typename graph_type, typename T>
1.41 + class preflow_push_max_flow {
1.42 +
1.43 + typedef typename graph_traits<graph_type>::node_iterator node_iterator;
1.44 + typedef typename graph_traits<graph_type>::each_node_iterator each_node_iterator;
1.45 + typedef typename graph_traits<graph_type>::out_edge_iterator out_edge_iterator;
1.46 + typedef typename graph_traits<graph_type>::in_edge_iterator in_edge_iterator;
1.47 +
1.48 + graph_type& G;
1.49 + node_iterator s;
1.50 + node_iterator t;
1.51 + edge_property_vector<graph_type, T>& capacity;
1.52 + T value;
1.53 + node_property_vector<graph_type, bool> mincutvector;
1.54 +
1.55 +
1.56 +
1.57 + public:
1.58 +
1.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) { }
1.60 +
1.61 +
1.62 + /*
1.63 + The run() function runs a modified version of the highest label preflow-push, which only
1.64 + finds a maximum preflow, hence giving the value of a maximum flow.
1.65 + */
1.66 + void run() {
1.67 +
1.68 + edge_property_vector<graph_type, T> flow(G, 0); //the flow value, 0 everywhere
1.69 + node_property_vector<graph_type, int> level(G); //level of node
1.70 + node_property_vector<graph_type, T> excess(G); //excess of node
1.71 +
1.72 + int n=number_of(G.first_node()); //number of nodes
1.73 + int b=n-2;
1.74 + /*b is a bound on the highest level of an active node. In the beginning it is at most n-2.*/
1.75 +
1.76 + std::vector<int> numb(n); //The number of nodes on level i < n.
1.77 +
1.78 + std::vector<std::stack<node_iterator> > stack(2*n-1); //Stack of the active nodes in level i.
1.79 +
1.80 +
1.81 +
1.82 + /*Reverse_bfs from t, to find the starting level.*/
1.83 +
1.84 + reverse_bfs<list_graph> bfs(G, t);
1.85 + bfs.run();
1.86 + for(each_node_iterator v=G.first_node(); v.valid(); ++v)
1.87 + {
1.88 + int dist=bfs.dist(v);
1.89 + level.put(v, dist);
1.90 + ++numb[dist];
1.91 + }
1.92 +
1.93 + /*The level of s is fixed to n*/
1.94 + level.put(s,n);
1.95 +
1.96 +
1.97 + /* Starting flow. It is everywhere 0 at the moment. */
1.98 +
1.99 + for(out_edge_iterator i=G.first_out_edge(s); i.valid(); ++i)
1.100 + {
1.101 + node_iterator w=G.head(i);
1.102 + flow.put(i, capacity.get(i));
1.103 + stack[bfs.dist(w)].push(w);
1.104 + excess.put(w, capacity.get(i));
1.105 + }
1.106 +
1.107 +
1.108 + /*
1.109 + End of preprocessing
1.110 + */
1.111 +
1.112 +
1.113 +
1.114 +
1.115 + /*
1.116 + Push/relabel on the highest level active nodes.
1.117 + */
1.118 +
1.119 + /*While there exists an active node.*/
1.120 + while (b) {
1.121 +
1.122 + /*We decrease the bound if there is no active node of level b.*/
1.123 + if (stack[b].empty()) {
1.124 + --b;
1.125 + } else {
1.126 +
1.127 + node_iterator w=stack[b].top(); //w is the highest label active node.
1.128 + stack[b].pop(); //We delete w from the stack.
1.129 +
1.130 + int newlevel=2*n-2; //In newlevel we maintain the next level of w.
1.131 +
1.132 + for(out_edge_iterator e=G.first_out_edge(w); e.valid(); ++e) {
1.133 + node_iterator v=G.head(e);
1.134 + /*e is the edge wv.*/
1.135 +
1.136 + if (flow.get(e)<capacity.get(e)) {
1.137 + /*e is an edge of the residual graph */
1.138 +
1.139 + if(level.get(w)==level.get(v)+1) {
1.140 + /*Push is allowed now*/
1.141 +
1.142 + if (capacity.get(e)-flow.get(e) > excess.get(w)) {
1.143 + /*A nonsaturating push.*/
1.144 +
1.145 + if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
1.146 + /*v becomes active.*/
1.147 +
1.148 + flow.put(e, flow.get(e)+excess.get(w));
1.149 + excess.put(v, excess.get(v)+excess.get(w));
1.150 + excess.put(w,0);
1.151 + //std::cout << w << " " << v <<" elore elen nonsat pump " << std::endl;
1.152 + break;
1.153 + } else {
1.154 + /*A saturating push.*/
1.155 +
1.156 + if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
1.157 + /*v becomes active.*/
1.158 +
1.159 + excess.put(v, excess.get(v)+capacity.get(e)-flow.get(e));
1.160 + excess.put(w, excess.get(w)-capacity.get(e)+flow.get(e));
1.161 + flow.put(e, capacity.get(e));
1.162 + //std::cout << w <<" " << v <<" elore elen sat pump " << std::endl;
1.163 + if (excess.get(w)==0) break;
1.164 + /*If w is not active any more, then we go on to the next node.*/
1.165 +
1.166 + } // if (capacity.get(e)-flow.get(e) > excess.get(w))
1.167 + } // if (level.get(w)==level.get(v)+1)
1.168 +
1.169 + else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);}
1.170 +
1.171 + } //if (flow.get(e)<capacity.get(e))
1.172 +
1.173 + } //for(out_edge_iterator e=G.first_out_edge(w); e.valid(); ++e)
1.174 +
1.175 +
1.176 +
1.177 + for(in_edge_iterator e=G.first_in_edge(w); e.valid(); ++e) {
1.178 + node_iterator v=G.tail(e);
1.179 + /*e is the edge vw.*/
1.180 +
1.181 + if (excess.get(w)==0) break;
1.182 + /*It may happen, that w became inactive in the first 'for' cycle.*/
1.183 +
1.184 + if(flow.get(e)>0) {
1.185 + /*e is an edge of the residual graph */
1.186 +
1.187 + if(level.get(w)==level.get(v)+1) {
1.188 + /*Push is allowed now*/
1.189 +
1.190 + if (flow.get(e) > excess.get(w)) {
1.191 + /*A nonsaturating push.*/
1.192 +
1.193 + if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
1.194 + /*v becomes active.*/
1.195 +
1.196 + flow.put(e, flow.get(e)-excess.get(w));
1.197 + excess.put(v, excess.get(v)+excess.get(w));
1.198 + excess.put(w,0);
1.199 + //std::cout << v << " " << w << " vissza elen nonsat pump " << std::endl;
1.200 + break;
1.201 + } else {
1.202 + /*A saturating push.*/
1.203 +
1.204 + if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
1.205 + /*v becomes active.*/
1.206 +
1.207 + flow.put(e,0);
1.208 + excess.put(v, excess.get(v)+flow.get(e));
1.209 + excess.put(w, excess.get(w)-flow.get(e));
1.210 + //std::cout << v <<" " << w << " vissza elen sat pump " << std::endl;
1.211 + if (excess.get(w)==0) { break;}
1.212 + } //if (flow.get(e) > excess.get(v))
1.213 + } //if(level.get(w)==level.get(v)+1)
1.214 +
1.215 + else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);}
1.216 + //std::cout << "Leveldecrease of node " << w << " to " << newlevel << std::endl;
1.217 +
1.218 + } //if (flow.get(e)>0)
1.219 +
1.220 + } //for in-edge
1.221 +
1.222 +
1.223 +
1.224 +
1.225 + /*
1.226 + Relabel
1.227 + */
1.228 + if (excess.get(w)>0) {
1.229 + /*Now newlevel <= n*/
1.230 +
1.231 + int l=level.get(w); //l is the old level of w.
1.232 + --numb[l];
1.233 +
1.234 + if (newlevel == n) {
1.235 + level.put(w,n);
1.236 +
1.237 + } else {
1.238 +
1.239 + if (numb[l]) {
1.240 + /*If the level of w remains nonempty.*/
1.241 +
1.242 + level.put(w,++newlevel);
1.243 + ++numb[newlevel];
1.244 + stack[newlevel].push(w);
1.245 + b=newlevel;
1.246 + } else {
1.247 + /*If the level of w gets empty.*/
1.248 +
1.249 + for (each_node_iterator v=G.first_node() ; v.valid() ; ++v) {
1.250 + if (level.get(v) >= l ) {
1.251 + level.put(v,n);
1.252 + }
1.253 + }
1.254 +
1.255 + for (int i=l+1 ; i!=n ; ++i) numb[i]=0;
1.256 + } //if (numb[l])
1.257 +
1.258 + } // if (newlevel = n)
1.259 +
1.260 + } // if (excess.get(w)>0)
1.261 +
1.262 +
1.263 + } //else
1.264 +
1.265 + } //while(b)
1.266 +
1.267 + value=excess.get(t);
1.268 + /*Max flow value.*/
1.269 +
1.270 +
1.271 +
1.272 + /*
1.273 + We find an empty level, e. The nodes above this level give
1.274 + a minimum cut.
1.275 + */
1.276 +
1.277 + int e=1;
1.278 +
1.279 + while(e) {
1.280 + if(numb[e]) ++e;
1.281 + else break;
1.282 + }
1.283 + for (each_node_iterator v=G.first_node(); v.valid(); ++v) {
1.284 + if (level.get(v) > e) mincutvector.put(v, true);
1.285 + }
1.286 +
1.287 +
1.288 + } // void run()
1.289 +
1.290 +
1.291 +
1.292 + /*
1.293 + Returns the maximum value of a flow.
1.294 + */
1.295 +
1.296 + T maxflow() {
1.297 + return value;
1.298 + }
1.299 +
1.300 +
1.301 +
1.302 + /*
1.303 + Returns a minimum cut.
1.304 + */
1.305 +
1.306 + node_property_vector<graph_type, bool> mincut() {
1.307 + return mincutvector;
1.308 + }
1.309 +
1.310 +
1.311 + };
1.312 +}//namespace marci
1.313 +#endif
1.314 +
1.315 +
1.316 +
1.317 +
1.318 +