1.1 --- a/src/work/jacint/preflow.h Thu Apr 29 15:01:52 2004 +0000
1.2 +++ b/src/work/jacint/preflow.h Thu Apr 29 15:58:34 2004 +0000
1.3 @@ -44,6 +44,11 @@
1.4 #include <stack>
1.5
1.6 #include <graph_wrapper.h>
1.7 +#include <bfs_iterator.h>
1.8 +#include <invalid.h>
1.9 +#include <maps.h>
1.10 +#include <for_each_macros.h>
1.11 +
1.12
1.13 namespace hugo {
1.14
1.15 @@ -111,9 +116,14 @@
1.16
1.17 bool augmentOnBlockingFlow2();
1.18
1.19 - //Returns the maximum value of a flow.
1.20 - Num flowValue() {
1.21 - return excess[t];
1.22 + /// Returns the actual flow value.
1.23 + /// More precisely, it returns the negative excess of s, thus
1.24 + /// this works also for preflows.
1.25 + Num flowValue() {
1.26 + Num a=0;
1.27 + FOR_EACH_INC_LOC(OutEdgeIt, e, *g, s) a+=(*flow)[e];
1.28 + FOR_EACH_INC_LOC(InEdgeIt, e, *g, s) a-=(*flow)[e];
1.29 + return a;
1.30 }
1.31
1.32 //should be used only between preflowPhase0 and preflowPhase1
1.33 @@ -433,7 +443,8 @@
1.34
1.35 void relabel(Node w, int newlevel, VecStack& active,
1.36 VecNode& level_list, NNMap& left,
1.37 - NNMap& right, int& b, int& k, bool what_heur ) {
1.38 + NNMap& right, int& b, int& k, bool what_heur )
1.39 + {
1.40
1.41 Num lev=level[w];
1.42
1.43 @@ -497,6 +508,28 @@
1.44 }
1.45
1.46 } //relabel
1.47 +
1.48 +
1.49 + template<typename MapGraphWrapper>
1.50 + class DistanceMap {
1.51 + protected:
1.52 + const MapGraphWrapper* g;
1.53 + typename MapGraphWrapper::template NodeMap<int> dist;
1.54 + public:
1.55 + DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
1.56 + void set(const typename MapGraphWrapper::Node& n, int a) {
1.57 + dist.set(n, a);
1.58 + }
1.59 + int operator[](const typename MapGraphWrapper::Node& n)
1.60 + { return dist[n]; }
1.61 +// int get(const typename MapGraphWrapper::Node& n) const {
1.62 +// return dist[n]; }
1.63 +// bool get(const typename MapGraphWrapper::Edge& e) const {
1.64 +// return (dist.get(g->tail(e))<dist.get(g->head(e))); }
1.65 + bool operator[](const typename MapGraphWrapper::Edge& e) const {
1.66 + return (dist[g->tail(e)]<dist[g->head(e)]);
1.67 + }
1.68 + };
1.69
1.70 };
1.71
1.72 @@ -674,8 +707,52 @@
1.73
1.74
1.75
1.76 + template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1.77 + bool Preflow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()
1.78 + {
1.79 + ResGW res_graph(*g, *capacity, *flow);
1.80 + bool _augment=false;
1.81 +
1.82 + //ReachedMap level(res_graph);
1.83 + FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
1.84 + BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
1.85 + bfs.pushAndSetReached(s);
1.86 +
1.87 + typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
1.88 + pred.set(s, INVALID);
1.89 +
1.90 + typename ResGW::template NodeMap<Num> free(res_graph);
1.91 +
1.92 + //searching for augmenting path
1.93 + while ( !bfs.finished() ) {
1.94 + ResGWOutEdgeIt e=bfs;
1.95 + if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
1.96 + Node v=res_graph.tail(e);
1.97 + Node w=res_graph.head(e);
1.98 + pred.set(w, e);
1.99 + if (res_graph.valid(pred[v])) {
1.100 + free.set(w, std::min(free[v], res_graph.resCap(e)));
1.101 + } else {
1.102 + free.set(w, res_graph.resCap(e));
1.103 + }
1.104 + if (res_graph.head(e)==t) { _augment=true; break; }
1.105 + }
1.106 +
1.107 + ++bfs;
1.108 + } //end of searching augmenting path
1.109
1.110 + if (_augment) {
1.111 + Node n=t;
1.112 + Num augment_value=free[t];
1.113 + while (res_graph.valid(pred[n])) {
1.114 + ResGWEdge e=pred[n];
1.115 + res_graph.augment(e, augment_value);
1.116 + n=res_graph.tail(e);
1.117 + }
1.118 + }
1.119
1.120 + return _augment;
1.121 + }
1.122
1.123
1.124
1.125 @@ -683,6 +760,253 @@
1.126
1.127
1.128
1.129 +
1.130 +
1.131 + template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1.132 + template<typename MutableGraph>
1.133 + bool Preflow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()
1.134 + {
1.135 + typedef MutableGraph MG;
1.136 + bool _augment=false;
1.137 +
1.138 + ResGW res_graph(*g, *capacity, *flow);
1.139 +
1.140 + //bfs for distances on the residual graph
1.141 + //ReachedMap level(res_graph);
1.142 + FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
1.143 + BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
1.144 + bfs.pushAndSetReached(s);
1.145 + typename ResGW::template NodeMap<int>
1.146 + dist(res_graph); //filled up with 0's
1.147 +
1.148 + //F will contain the physical copy of the residual graph
1.149 + //with the set of edges which are on shortest paths
1.150 + MG F;
1.151 + typename ResGW::template NodeMap<typename MG::Node>
1.152 + res_graph_to_F(res_graph);
1.153 + {
1.154 + typename ResGW::NodeIt n;
1.155 + for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
1.156 + res_graph_to_F.set(n, F.addNode());
1.157 + }
1.158 + }
1.159 +
1.160 + typename MG::Node sF=res_graph_to_F[s];
1.161 + typename MG::Node tF=res_graph_to_F[t];
1.162 + typename MG::template EdgeMap<ResGWEdge> original_edge(F);
1.163 + typename MG::template EdgeMap<Num> residual_capacity(F);
1.164 +
1.165 + while ( !bfs.finished() ) {
1.166 + ResGWOutEdgeIt e=bfs;
1.167 + if (res_graph.valid(e)) {
1.168 + if (bfs.isBNodeNewlyReached()) {
1.169 + dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
1.170 + typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
1.171 + original_edge.update();
1.172 + original_edge.set(f, e);
1.173 + residual_capacity.update();
1.174 + residual_capacity.set(f, res_graph.resCap(e));
1.175 + } else {
1.176 + if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
1.177 + typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
1.178 + original_edge.update();
1.179 + original_edge.set(f, e);
1.180 + residual_capacity.update();
1.181 + residual_capacity.set(f, res_graph.resCap(e));
1.182 + }
1.183 + }
1.184 + }
1.185 + ++bfs;
1.186 + } //computing distances from s in the residual graph
1.187 +
1.188 + bool __augment=true;
1.189 +
1.190 + while (__augment) {
1.191 + __augment=false;
1.192 + //computing blocking flow with dfs
1.193 + DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
1.194 + typename MG::template NodeMap<typename MG::Edge> pred(F);
1.195 + pred.set(sF, INVALID);
1.196 + //invalid iterators for sources
1.197 +
1.198 + typename MG::template NodeMap<Num> free(F);
1.199 +
1.200 + dfs.pushAndSetReached(sF);
1.201 + while (!dfs.finished()) {
1.202 + ++dfs;
1.203 + if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
1.204 + if (dfs.isBNodeNewlyReached()) {
1.205 + typename MG::Node v=F.aNode(dfs);
1.206 + typename MG::Node w=F.bNode(dfs);
1.207 + pred.set(w, dfs);
1.208 + if (F.valid(pred[v])) {
1.209 + free.set(w, std::min(free[v], residual_capacity[dfs]));
1.210 + } else {
1.211 + free.set(w, residual_capacity[dfs]);
1.212 + }
1.213 + if (w==tF) {
1.214 + __augment=true;
1.215 + _augment=true;
1.216 + break;
1.217 + }
1.218 +
1.219 + } else {
1.220 + F.erase(/*typename MG::OutEdgeIt*/(dfs));
1.221 + }
1.222 + }
1.223 + }
1.224 +
1.225 + if (__augment) {
1.226 + typename MG::Node n=tF;
1.227 + Num augment_value=free[tF];
1.228 + while (F.valid(pred[n])) {
1.229 + typename MG::Edge e=pred[n];
1.230 + res_graph.augment(original_edge[e], augment_value);
1.231 + n=F.tail(e);
1.232 + if (residual_capacity[e]==augment_value)
1.233 + F.erase(e);
1.234 + else
1.235 + residual_capacity.set(e, residual_capacity[e]-augment_value);
1.236 + }
1.237 + }
1.238 +
1.239 + }
1.240 +
1.241 + return _augment;
1.242 + }
1.243 +
1.244 +
1.245 +
1.246 +
1.247 +
1.248 +
1.249 + template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1.250 + bool Preflow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()
1.251 + {
1.252 + bool _augment=false;
1.253 +
1.254 + ResGW res_graph(*g, *capacity, *flow);
1.255 +
1.256 + //ReachedMap level(res_graph);
1.257 + FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
1.258 + BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
1.259 +
1.260 + bfs.pushAndSetReached(s);
1.261 + DistanceMap<ResGW> dist(res_graph);
1.262 + while ( !bfs.finished() ) {
1.263 + ResGWOutEdgeIt e=bfs;
1.264 + if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
1.265 + dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
1.266 + }
1.267 + ++bfs;
1.268 + } //computing distances from s in the residual graph
1.269 +
1.270 + //Subgraph containing the edges on some shortest paths
1.271 + ConstMap<typename ResGW::Node, bool> true_map(true);
1.272 + typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
1.273 + DistanceMap<ResGW> > FilterResGW;
1.274 + FilterResGW filter_res_graph(res_graph, true_map, dist);
1.275 +
1.276 + //Subgraph, which is able to delete edges which are already
1.277 + //met by the dfs
1.278 + typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt>
1.279 + first_out_edges(filter_res_graph);
1.280 + typename FilterResGW::NodeIt v;
1.281 + for(filter_res_graph.first(v); filter_res_graph.valid(v);
1.282 + filter_res_graph.next(v))
1.283 + {
1.284 + typename FilterResGW::OutEdgeIt e;
1.285 + filter_res_graph.first(e, v);
1.286 + first_out_edges.set(v, e);
1.287 + }
1.288 + typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
1.289 + template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
1.290 + ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
1.291 +
1.292 + bool __augment=true;
1.293 +
1.294 + while (__augment) {
1.295 +
1.296 + __augment=false;
1.297 + //computing blocking flow with dfs
1.298 + DfsIterator< ErasingResGW,
1.299 + typename ErasingResGW::template NodeMap<bool> >
1.300 + dfs(erasing_res_graph);
1.301 + typename ErasingResGW::
1.302 + template NodeMap<typename ErasingResGW::OutEdgeIt>
1.303 + pred(erasing_res_graph);
1.304 + pred.set(s, INVALID);
1.305 + //invalid iterators for sources
1.306 +
1.307 + typename ErasingResGW::template NodeMap<Num>
1.308 + free1(erasing_res_graph);
1.309 +
1.310 + dfs.pushAndSetReached(
1.311 + typename ErasingResGW::Node(
1.312 + typename FilterResGW::Node(
1.313 + typename ResGW::Node(s)
1.314 + )
1.315 + )
1.316 + );
1.317 + while (!dfs.finished()) {
1.318 + ++dfs;
1.319 + if (erasing_res_graph.valid(
1.320 + typename ErasingResGW::OutEdgeIt(dfs)))
1.321 + {
1.322 + if (dfs.isBNodeNewlyReached()) {
1.323 +
1.324 + typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
1.325 + typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
1.326 +
1.327 + pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
1.328 + if (erasing_res_graph.valid(pred[v])) {
1.329 + free1.set(w, std::min(free1[v], res_graph.resCap(
1.330 + typename ErasingResGW::OutEdgeIt(dfs))));
1.331 + } else {
1.332 + free1.set(w, res_graph.resCap(
1.333 + typename ErasingResGW::OutEdgeIt(dfs)));
1.334 + }
1.335 +
1.336 + if (w==t) {
1.337 + __augment=true;
1.338 + _augment=true;
1.339 + break;
1.340 + }
1.341 + } else {
1.342 + erasing_res_graph.erase(dfs);
1.343 + }
1.344 + }
1.345 + }
1.346 +
1.347 + if (__augment) {
1.348 + typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t));
1.349 +// typename ResGW::NodeMap<Num> a(res_graph);
1.350 +// typename ResGW::Node b;
1.351 +// Num j=a[b];
1.352 +// typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
1.353 +// typename FilterResGW::Node b1;
1.354 +// Num j1=a1[b1];
1.355 +// typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
1.356 +// typename ErasingResGW::Node b2;
1.357 +// Num j2=a2[b2];
1.358 + Num augment_value=free1[n];
1.359 + while (erasing_res_graph.valid(pred[n])) {
1.360 + typename ErasingResGW::OutEdgeIt e=pred[n];
1.361 + res_graph.augment(e, augment_value);
1.362 + n=erasing_res_graph.tail(e);
1.363 + if (res_graph.resCap(e)==0)
1.364 + erasing_res_graph.erase(e);
1.365 + }
1.366 + }
1.367 +
1.368 + } //while (__augment)
1.369 +
1.370 + return _augment;
1.371 + }
1.372 +
1.373 +
1.374 +
1.375 +
1.376 } //namespace hugo
1.377
1.378 #endif //HUGO_PREFLOW_H
2.1 --- a/src/work/marci/edmonds_karp_demo.cc Thu Apr 29 15:01:52 2004 +0000
2.2 +++ b/src/work/marci/edmonds_karp_demo.cc Thu Apr 29 15:58:34 2004 +0000
2.3 @@ -5,11 +5,11 @@
2.4 #include <list_graph.h>
2.5 #include <smart_graph.h>
2.6 #include <dimacs.h>
2.7 -#include <edmonds_karp.h>
2.8 +//#include <edmonds_karp.h>
2.9 #include <time_measure.h>
2.10 //#include <graph_wrapper.h>
2.11 #include <preflow.h>
2.12 -#include <preflow_res.h>
2.13 +//#include <preflow_res.h>
2.14 #include <for_each_macros.h>
2.15
2.16 using namespace hugo;
2.17 @@ -72,12 +72,12 @@
2.18 Graph::EdgeMap<int> flow(G); //0 flow
2.19 Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
2.20 pre_flow_test(G, s, t, cap, flow/*, true*/);
2.21 - Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
2.22 - pre_flow_ize(G, s, t, cap, flow/*, false*/);
2.23 + // Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
2.24 + // pre_flow_ize(G, s, t, cap, flow/*, false*/);
2.25 // PreflowRes<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
2.26 // pre_flow_res(G, s, t, cap, flow/*, true*/);
2.27 - MaxFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
2.28 - max_flow_test(G, s, t, cap, flow);
2.29 + //MaxFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
2.30 + // max_flow_test(G, s, t, cap, flow);
2.31
2.32 {
2.33 std::cout << "preflow ..." << std::endl;
2.34 @@ -91,9 +91,9 @@
2.35 std::cout << "preflow ..." << std::endl;
2.36 FOR_EACH_LOC(Graph::EdgeIt, e, G) flow.set(e, 0);
2.37 ts.reset();
2.38 - pre_flow_ize.preflow(Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
2.39 + pre_flow_test.preflow(Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
2.40 std::cout << "elapsed time: " << ts << std::endl;
2.41 - std::cout << "flow value: "<< pre_flow_ize.flowValue() << std::endl;
2.42 + std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
2.43 }
2.44
2.45 // {
2.46 @@ -110,32 +110,32 @@
2.47 FOR_EACH_LOC(Graph::EdgeIt, e, G) flow.set(e, 0);
2.48 ts.reset();
2.49 int i=0;
2.50 - while (max_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
2.51 + while (pre_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
2.52 std::cout << "elapsed time: " << ts << std::endl;
2.53 std::cout << "number of augmentation phases: " << i << std::endl;
2.54 - std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
2.55 + std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
2.56 }
2.57
2.58 - {
2.59 - std::cout << "faster physical blocking flow augmentation ..." << std::endl;
2.60 - FOR_EACH_LOC(Graph::EdgeIt, e, G) flow.set(e, 0);
2.61 - ts.reset();
2.62 - int i=0;
2.63 - while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
2.64 - std::cout << "elapsed time: " << ts << std::endl;
2.65 - std::cout << "number of augmentation phases: " << i << std::endl;
2.66 - std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
2.67 - }
2.68 +// {
2.69 +// std::cout << "faster physical blocking flow augmentation ..." << std::endl;
2.70 +// FOR_EACH_LOC(Graph::EdgeIt, e, G) flow.set(e, 0);
2.71 +// ts.reset();
2.72 +// int i=0;
2.73 +// while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
2.74 +// std::cout << "elapsed time: " << ts << std::endl;
2.75 +// std::cout << "number of augmentation phases: " << i << std::endl;
2.76 +// std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
2.77 +// }
2.78
2.79 {
2.80 std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
2.81 FOR_EACH_LOC(Graph::EdgeIt, e, G) flow.set(e, 0);
2.82 ts.reset();
2.83 int i=0;
2.84 - while (max_flow_test.augmentOnBlockingFlow2()) { ++i; }
2.85 + while (pre_flow_test.augmentOnBlockingFlow2()) { ++i; }
2.86 std::cout << "elapsed time: " << ts << std::endl;
2.87 std::cout << "number of augmentation phases: " << i << std::endl;
2.88 - std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
2.89 + std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
2.90 }
2.91
2.92 {
2.93 @@ -143,10 +143,10 @@
2.94 FOR_EACH_LOC(Graph::EdgeIt, e, G) flow.set(e, 0);
2.95 ts.reset();
2.96 int i=0;
2.97 - while (max_flow_test.augmentOnShortestPath()) { ++i; }
2.98 + while (pre_flow_test.augmentOnShortestPath()) { ++i; }
2.99 std::cout << "elapsed time: " << ts << std::endl;
2.100 std::cout << "number of augmentation phases: " << i << std::endl;
2.101 - std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
2.102 + std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
2.103 }
2.104
2.105